CC BY
120
_МЕЖДУНАРОДНЫЙ НАУЧНЫЙ ЖУРНАЛ «СИМВОЛ НАУКИ» №10/2015 ISSN 2410-700Х_
УДК 518:517.949.8; 533.6.011
Ворожцов Евгений Васильевич
вед. науч. тотрудник ИТПМ СО РАН, г. Новосибирск, РФ E-mail: [email protected] Шапеев Василий Павлович главн. науч. сотрудник ИТПМ СО РАН, г. Новосибирск, РФ E-mail: [email protected]
ОБ УСКОРЕНИИ ИТЕРАЦИОННЫХ ПРОЦЕССОВ РЕШЕНИЯ КРАЕВЫХ ЗАДАЧ КОМБИНИРОВАНИЕМ МЕТОДОВ КРЫЛОВА И ФЕДОРЕНКО1
Аннотация
Реализованы несколько новых вариантов известных способов ускорения процессов итерационного решения дискретных задач, которые возникают при решении численными методами краевых задач для уравнений с частными производными (PDE). В частности, предложено и реализовано комбинированное применение операции продолжения метода Федоренко и метода Крылова. В качестве объекта численных экспериментов взята дискретная задача, к которой сводится решение уравнений Навье-Стокса методом коллокаций и наименьших невязок (КНН). Показано значительное (в сотни раз) ускорение решения задачи в случае комбинирования указанных двух способов. Обсуждается, какие специфические свойства метода КНН повлияли на результат. Выводы подкреплены результатами многих численных экспериментов.
Ключевые слова
ускорение итерационных процессов, подпространства Крылова, многосеточные алгоритмы, метод коллокаций и наименьших невязок, уравнения Навье-Стокса.
1. Введение
В большинстве численных методов решения краевых задач для нелинейных PDE используются линеаризация и итерации по нелинейности. Путем того или иного способа дискретизации — перехода от исходной задачи для PDE (дифференциальной задачи) к ее дискретной аппроксимации (дискретной задаче) приближенное решение дифференциальной задачи в большинстве методов сводится к решению систем линейных алгебраических уравнений (СЛАУ). При стремлении достичь хорошей точности численного решения вводится большое количество неизвестных дискретной задачи. Так, например, зачастую в численных методах для этого вводится большое количество подлежащих вычислению неизвестных значений сеточных функций на сетках с большим количеством узлов. Это приводит к необходимости решать большого размера СЛАУ. Требуемое для этого количество арифметических операций зачастую приводит к большому времени решения дискретной задачи на ЭВМ и/или к накоплению большой ошибки округлений. Каждое из этих двух обстоятельств может сделать неприемлемым применение прямых методов для решения конкретной дискретной задачи. По этой причине востребованы итерационные методы решения СЛАУ, поскольку прямые методы их решения зачастую требуют существенно больше арифметических операций при однократном применении, чем итерационные методы в течение одной итерации. Если итерации сходятся, их можно многократно повторить до достижения некоторой точности решения задачи. При этом проблема сокращения времени решения задачи остается.
Поиск новых численных алгоритмов и более эффективных реализаций на ЭВМ новых и существующих методов может расширить возможности персональных компьютеров, а также существенно усилить важное направление ускорения решения задач с применением распараллеливания алгоритмов.
1 Работа выполнена при частичной поддержке РФФИ (грант 13-01-00277).
_МЕЖДУНАРОДНЫЙ НАУЧНЫЙ ЖУРНАЛ «СИМВОЛ НАУКИ» №10/2015 ISSN 2410-700Х_
В численных методах широкое применение нашли следующие подходы к ускорению итерационных процессов решения задач для PDE: ускорение сходимости решения СЛАУ с применением подпространств Крылова, ускорение сходимости с помощью многосеточных алгоритмов, ускорение сходимости с помощью предобусловливателей (в данной работе не рассматривается).
Ускорение сходимости с помощью подпространств Крылова уже давно применяется при численном решении для PDE. В частности, в работе [1] описан итерационный метод Ньютона-Крылова для ускорения итераций численного решения стационарных уравнений Навье-Стокса. В [2] дан обзор методов Ньютона— Крылова. В [3] возникающая в результате дискретизации исходных PDE-уравнений СЛАУ решалась с применением гибкого алгоритма GMRES, являющегося вариантом реализации метода Крылова. Различные его варианты приведены в [4].
Примененный здесь способ ускорения итерационных процессов с использованием подпространств Крылова является новым вариантом известного и описанного ранее в [5, 6]. Он также применялся ранее для ускорения сходимости метода КНН решения двумерных уравнений Навье-Стокса [7].
Другой способ ускорения сходимости итерационного процесса решения стационарных уравнений Навье-Стокса, который применяется в настоящей статье, основан на многосеточных алгоритмах. Основная идея этих алгоритмов была предложена Р.П. Федоренко [8]. В настоящее время многосеточные алгоритмы известны как очень эффективные методы численного решения разнообразных задач математической физики. Их основная идея состоит в том, чтобы ускорять сходимость основного итерационного процесса путем решения задачи последовательно на сетках различного разрешения [9].
В данной работе рассматривается комбинированное применение методов Крылова и Федоренко для ускорения итерационного решения краевых задач. Численные эксперименты для наглядности проведены при решении уравнений Навье-Стокса. Известно, что вне зависимости от применяемых численных методов их решения при больших числах Рейнольдса (> 1000) для достижения хорошей точности возникает необходимость решать большие СЛАУ и применять итерационные методы, требующие большие временные затраты на ЭВМ. Здесь, в отличие от [10], для простоты взяты двумерные уравнения.
Следует отметить, что комбинирование методов Крылова и Федоренко ранее применялось в ряде численных методов, существенно отличающихся от метода КНН. Оно применялось в случаях, когда СЛАУ возникала в результате дискретизации уравнений Навье-Стокса с помощью конечно-разностных методов и методов конечного объема [3, 11-15]. В частности, техника ускорения итерационных процессов путем последовательного включения многосеточного и GMRES алгоритмов была реализована в [14]. В результате было достигнуто ускорение до 25 раз по сравнению со стандартным нелинейным многосеточным алгоритмом. В работах [16-19] описаны некоторые варианты комбинирования методов Крылова и Федоренко в рамках различных вариантов метода конечных элементов (МКЭ).
Метод КНН обладает рядом достоинств. Ввиду наличия кусочно-аналитического решения он удобен для применения в областях с криволинейными границами и на адаптивных сетках [20]; в нем также относительно легко формулируются варианты схем повышенного порядка точности, за счет увеличения степени базисных полиномов [21]; метод удобен для распараллеливания благодаря тому, что в нем реализован альтернирующий метод Шварца. Минимизация функционала невязки уравнений дискретной задачи и условия гладкого сопряжения кусков аналитического ее решения на границах между соседними ячейками позволяют в методе КНН подавлять не соответствующие решению дифференциальной задачи "посторонние" осцилляции численного решения [22]. В нем есть возможность до решения дискретной задачи за счет выбора коэффициентов разложения решения в базисе функционального пространства удовлетворить (с точностью выбранного порядка аппроксимации) тождественно часть линейных уравнений задачи. Например, в использованном здесь варианте метода это сделано относительно уравнения неразрывности в системе уравнений Навье-Стокса. Решение его совместно с уравнениями импульса представляет собой определенные трудности для отдельных численных методов.
Численные эксперименты по ускорению сходимости итераций проведены здесь при решении уравнений Навье-Стокса методом КНН — вариантом метода коллокаций и наименьших квадратов (КНК), который появился относительно недавно в результате поиска новых методов решения PDE [23-25]. Этому
_МЕЖДУНАРОДНЫЙ НАУЧНЫЙ ЖУРНАЛ «СИМВОЛ НАУКИ» №10/2015 ISSN 2410-700Х_
предшествовал ряд работ, в которых методом коллокаций эффективно решались краевые задачи для обыкновенных дифференциальных уравнений (ОДУ). Достаточно простой в реализации метод коллокаций, который сравнительно легко может применяться на нерегулярных сетках, дает хорошие результаты при решении краевых задач для обыкновенных дифференциальных уравнений и хорошо обусловленных задач для эллиптических и параболических уравнений. Но попытки применить его к решению задач с особенностями, например, в виде узких пограничных слоев с большими градиентами решения приводят к плохо обусловленным СЛАУ в соответствующей дискретной задаче. Показательными в этом отношении являются краевые задачи для уравнений Навье-Стокса при больших числах Рейнольдса, которые представляют аналогичные трудности и для некоторых других численных методов. В методе КНК метод коллокаций сочетается с открытым Лежандром и независимо Гауссом методом наименьших квадратов (МНК). Последний применяется в решении задачи построения аппроксиманта функции (ПАФ), заданной дискретно с некоторой погрешностью в значениях. В методах КНК и КНН для построения приближенного решения исходной дифференциальной задачи, как и в задаче ПАФ, сначала выписывается переопределенная СЛАУ. Затем с помощью МНК отыскивается некоторое ее "псевдорешение", которое является приближенным численным решением исходной задачи. Вследствие применения МНК на соответствующем псевдорешении переопределенной СЛАУ достигается минимум функционала невязки, составленного из невязок всех ее уравнений. Последнее опосредованно связано с погрешностью решения исходной задачи. Вследствие минимизации функционала невязки в дискретной задаче, как показывают численные эксперименты, методами КНК и КНН можно решать задачи с особенностями в решении, для которых неприменим метод коллокаций. Показателен в этом отношении метод КНН, в котором для решения переопределенных СЛАУ используется ортогональный метод решения переопределенной СЛАУ вместо МНК, применяемого в методе КНК. Это связано с тем, что в первом случае в сравнении со вторым не ухудшается обусловленность матрицы СЛАУ дискретной задачи [10].
В работе [7] метод КНН был обобщен на случай нерегулярной сетки с прямоугольными ячейками в плоскости двух пространственных переменных. В [7] представлены варианты метода КНН, имеющие порядок точности от 2 до 8. В [7] возможности метода КНН были показаны на решении известной эталонной двумерной задачи о динамике вязкой жидкости в каверне с движущейся крышкой. Все возрастающей точности результаты по ее решению вот уже на протяжении более, чем тридцати лет, публикуют исследователи. В работах [7, 26] было проведено детальное сравнение ее численных решений и показано, что полученные методом КНН результаты не уступают по точности наиболее точным результатам, полученным другими авторами (см. [21, 27]) с помощью различных численных методов высокой точности.
В [10, 28] метод КНН применен для решения трехмерных уравнений Навье--Стокса несжимаемой жидкости на пространственной расчетной сетке из кубических ячеек. В качестве примера приведены результаты расчетов течения в кубической каверне с движущейся крышкой при Re = 100 и Re =1000, показано их сравнение с известными наиболее точными расчетами.
Комбинирование в методе КНН метода Крылова с операцией продолжения на многосеточном комплексе, являющемся составной частью метода Федоренко, впервые реализовано в данной работе. Следовало ожидать, что такая комбинация даст, как и в перечисленных выше случаях, существенное ускорение решения задачи на ЭВМ. Исследование эффективности применения комбинации методов Федоренко и нового варианта алгоритма Крылова для ускорения решения PDE составляет содержание представленной работы.
1. Описание варианта метода КНН численного решения двумерных уравнений Навье--Стокса
1.1. Постановка дифференциальной задачи
Рассмотрим краевую задачу для системы стационарных уравнений Навье-Стокса
(V -V)V + Vp = — SV - f, (1)
Re
divV = 0, (x, x2) eQ, (2)
По = g, (3)
которые описывают течение вязкой нетеплопроводной несжимаемой жидкости в области О с границей дО. В уравнениях(1) и (2) xi, Х2 — декартовы пространственные координаты, V = V(xi, x2) — вектор скорости, имеющий компоненты vl (Xl, X2),V2 (Xl, X2) вдоль осей xi, X2, соответственно; p = p(xi,x2)
д2 d2
— давление, f = f1, /2) — заданная вектор-функция, Re — число Рейнольдса, А —
дх1 дх 2
У'
/тл ТТЛ- ^ д
(к • V) = V--+ V-• Система из трех уравнений (1) и (2) решается при краевых условиях Дирихле
дх^ дх^
(3), где § = § (Х1, Х2 ) = (§, §2 ) — заданная вектор-функция. На давление накладывается условие
Jdx]dx2 = 0, (4)
которое справедливо при отсутствии источников и стоков в области О [7]. 1.2. Локальные координаты и базисные функции
В дальнейшем в качестве области решения задачи берется О квадрат
О = {(х,х2),0 < X < X, 1 = 1,2}, (5)
где X >0 — заданная длина стороны квадрата. В данной задаче (1)-(5) расчетная область (5) дискретизируется сеткой с квадратными ячейками С2/ .. /, / = 1,...,/ . где / — заданное число ячеек на
каждой из осей х\, Х2, I > 1 . Решение ищется на этой сетке в виде кусочно-гладкой функции. Для записи формул метода КНН удобно ввести локальные координаты у, у в каждой ячейке О,, . Зависимость
локальных координат от глобальных пространственных переменных Х1, Х2 задается соотношениями Ут = (Хт — Хт1 .)/ И, т = 1, 2, где Хт,ц — значение координаты хт в геометрическом центре ячейки О, ., а к — половина длины стороны квадратной ячейки О . Следовательно, локальные координаты изменяются в интервале ут е[—1,1] . Введем обозначения и(у, у ) = (щ, Щ ) = К(ИУ1 + Хц, j, Иу2 + х2.^) , р( У1, У2) = р(Иу + , Иу2 + х24и). В результате этой замены переменных уравнения Навье-Стокса принимают следующий вид:
Аи - Reh
Г дит дит др Л щ —- + щ —- +
V дУ д^2 дУт У
1 h
^дщ дм? ^ 1 t 2
дУ2 у
—Re- h2fm т — 12; (6)
— 0, (7)
Л - д2 д2
где А —-— +--—. Линеаризуем по Ньютону уравнения Навье-Стокса (6)
дУ дУ2
/+1 - (Re- h) (иЩ^1 + ufus + usus+1 + ufus
т V / \ 1 т, y 1 т, y 2 т, yn 2 т,
ДЩ+1 — (Де И)Щи^1 + Щ+Щ + и!и5+1 + щ5+1и5 + р5+1) = ^ , т = 1,2, (8)
т V \ 1 т,у 1 т,у 2 т,у2 2 т,у2 Ут / т' •> •> \ >
где 5 — номер итерации по нелинейности, = 0,1,2,..., И^,, л — известное приближение к решению на 5-ой итерации, начиная от выбранного на начального приближения с индексом 5 = 0,
К —Re
т
h2 fm - h (u;Vv + uU v |l, м , — дм / сУ,, pv — др / дк,, l, m = 1, 2.
•J т 11 т,y1 2 т,y^ I P т,y^ т s i ^ г y^ г J т' ' '
Приближенное решение в каждой ячейке а, ищется в виде линейной комбинации базисных векторов-функций <:
(<, ,р )т = ХЬ }<, (9)
/
где верхний индекс Т обозначает операцию транспонирования. В данном варианте метода решение задачи проектируется в пространство многочленов. Если число ячеек в области а больше единицы, то решение дискретной задачи является кусочно полиномиальным, т.е. состоит из кусков, связанных с началом локальной системы координат в каждой ячейке. В данной работе для аппроксимации составляющих скорости использованы многочлены второй степени по переменным у, у2, а для аппроксимации давления - многочлены первой степени.
Таблица 1
Вид базисных функций <
/ 1 2 3 4 5 6 7 8 9 10 11 12
1 У1 У2 у2 -2 У1У2 У2 0 0 0 0 0 0
0 ~У2 0 -2 у У2 у2 0 1 У1 у2 0 0 0
0 0 0 0 0 0 0 0 0 1 У1 У2
Всего базисных функций в выбранном пространстве пятнадцать. Поскольку в простом по записи уравнении неразрывности коэффициенты константы, то его легко удовлетворить за счет выбора базисных полиномов <. Нетрудно установить, что для этого требуется удовлетворить ими три линейных
соотношения. В итоге из первоначальных пятнадцати базисных полиномов независимыми останутся только двенадцать. Они приведены в табл. 1. Их совокупность можно назвать соленоидальным базисом, так как =0. В таком базисе уравнение неразрывности решением дискретной задачи с точностью до
выбранного порядка аппроксимации удовлетворяется тождественно в каждой ячейке. Уменьшение числа неизвестных коэффициентов в искомом решении (9) приводит к заметной экономии машинного времени при решении задачи.
1.3. Уравнения дискретной задачи
Метод КНН существенно отличается от других методов решения краевых задач для дифференциальных уравнений тем, что решение дискретной задачи определяется решением переопределенной СЛАУ. При этом требуется, чтобы на псевдорешении определяющей СЛАУ был минимум функционала невязки ее уравнений. Для записи уравнений дискретной задачи, определяющих
решение в ячейке а, , зададим в ней точки коллокации. Количество этих точек и их расположение внутри ячейки может варьироваться в различных вариантах метода. В данной работе были реализованы четыре варианта задания координат точек коллокации. Обозначим через Nчисло точек коллокации внутри каждой ячейки. В случае, когда N = 2, координаты точек коллокации таковы: (с,с), (—С,С), где С
— задаваемое пользователем значение в интервале 0 < С < 1. На рис. 1 точки коллокации показаны черными кружками, при этом С = 1 / 2. На рис. 1, (б) показано расположение точек коллокации при N = 4, при этом локальные координаты этих точек имеют вид (¿С, ¿с) , и при С = 1 / 2 указанные точки лежат в серединах полудиагоналей квадратной ячейки. В случае N = 8 (рис. 1, (в)) координаты точек коллокации задавались следующим образом: расположение первых четырех точек было взято таким же, как на рис. 1, (б), а координаты следующих четырех точек задавались по формулам (+С,0) , (0, +с) .
(а)
(б)
(в)
Рисунок 1 - Варианты задания точек коллокации и согласования: (а) Nc — 2, Nm = 4; (б) Nc = 4, Nm = 8
(в) Nc = 8.
Подставляя (9) в (8) и выписывая полученные соотношения в каждой точке коллокации с численными значениями ее координат, получим 2Nc уравнений коллокаций дискретной задачи:
12
v = l ,...,2NC. (Ю)
m—1
Они являются линейными алгебраическими уравнениями относительно искомых величин Ъ/ . , в (9).
Дополним систему уравнений дискретной задачи в ячейке а условиями согласования с решениями
дискретной задачи, взятыми во всех соседних с ней ячейках. Эти условия, которые не должны противоречить определению решения дифференциальной задачи, запишем в отдельных точках (точках
согласования) на сторонах ячейки а, общих с ее соседними ячейками. В данном случае условия согласования — это условия гладкого сопряжения решения дискретной задачи в точках согласования:
, д(и+)п . +.и . д(и— )п . ,и . д(и+) . +.г . д(и—) . к——— + (и ) = к——— + (и ) , к——— + (и ) = к——— + (и У, (11)
дп
дп
1 д , Здесь h — — h
дп
дд п--ь п
V
1 дхх
дх
— п
д
п
дп p+ — p~ •
д
дп
(12)
2 У
1 ду1 2 ду.
, п — (п, П ) — внешняя нормаль к стороне
ячейки а j, (")", (")Г — нормальные и касательные составляющие вектора скорости по отношению к + —
стороне ячейки, и , и — пределы функции и при стремлении ее аргументов к точке согласования изнутри и снаружи ячейки а .
Для единственности определения давления в решении задаем его значение в одной точке области либо аппроксимируем условие (4) по формуле
\ 1 { \
. (13)
7 У к V 7 У
Т*
Здесь I — интеграл по всей области, рассчитанный как сумма интегралов по каждой ячейке на
*
предыдущей итерации, р — давление в ячейке с предыдущей итерации.
1 ff pdQn 1 —1 (-I * + f p*dQii
h lJV 1J J h { 4/ 1J,
_МЕЖДУНАРОДНЫЙ НАУЧНЫЙ ЖУРНАЛ «СИМВОЛ НАУКИ» №10/2015 ISSN 2410-700Х_
Отметим, что в условиях (11) в записи систем уравнений, определяющих куски решения в двух соседних ячейках, направления внешних нормалей к этим ячейкам на их общей стороне противоположны. Поэтому в формулах (11) для одной из этих ячеек первые слагаемые войдут с одним знаком, а в соседней ячейке с противоположным. Благодаря этому, уравнения (11) не повторяются в совокупной СЛАУ, определяющей численное решение глобальной задачи и состоящей из всех уравнений, выписанных для каждой ячейки.
Точки, в которых записываются уравнения (11), называются точками согласования. Так же, как и в случае точек коллокации, количество точек согласования решения и их расположение на каждой стороне ячейки может быть различным в различных конкретных вариантах метода. Обозначим через Nm число точек согласования для составляющих вектора скорости на сторонах каждой ячейки. Некоторые из способов задания точек согласования, которые были реализованы в данной работе, показаны на рис. 0 малыми квадратами. На рис. 1, (а) координаты четырех точек согласования для составляющих вектора скорости заданы формулами (+1,0) , (0,¿1). В случае, изображенном на рис. 1, (б) 8 точек согласования
с координатами (+1, — £), (+1,£), (—£, + 1), (£, + 1), где параметр £ задается вычислителем в
интервале 0 < £ < 1. В частности, на рис. 0 расположение точек согласования показано для случая, когда
£ — 1/2.
Используя представление решения (9), подставим координаты этих точек в каждое из трех условий согласования (11). Из первых двух условий получим 2Nm линейных алгебраических уравнений для компонент скорости.
Условия согласования для давления (12) задаются в четырех точках с координатами (+1,0), (0, + 1). Эти точки показаны на рис. 1, (а). После подстановки представления (9) в (12) из этих условий
также получаются четыре линейных алгебраических уравнения (согласования).
Чтобы обеспечить единственность решения для давления, его значение задавалось в начале
глобальной системы координат (в данном работе в вершине ячейки О, ) или использовалось условие (13).
В тестовых решениях бралось точное значение давления в этой точке. Таким образом, условия
согласования составляющих скорости и давления дают в общей сложности 2Nт + 4 + линейных
относительно искомых b j i алгебраических уравнений в каждой ячейке (i, j), где S^ - символ
Кронекера, S^ — 1 при i — j и SS — 0 при i Ф j .
Если сторона ячейки совпадает с границей области О, то в соответствующих точках вместо условий согласования для решения дискретной задачи выписываются граничные условия: u — gm, т — 1,2.
Объединяя уравнения коллокаций, согласования и уравнения, полученные из краевых условий, если ячейка О, граничная, в каждой ячейке получим СЛАУ вида
.с1, (14)
где Ь^ = (/^'i• •^//'п )7 • Решение системы (14) определяет локально в данной ячейке Qг- . решение глобальной дискретной задачи как ее кусок, согласованный с соседними кусками. Матрица A j
содержит 2Nc + 2Nm + 4 + S-Sj строк и 12 столбцов. В вариантах, исследованных в данной работе,
система (14) переопределена. Для вывода уравнений этой СЛАУ полезно применение систем компьютерной алгебры. Здесь использовалась система Mathematica.
Здесь для решения СЛАУ дискретной задачи применялся процесс, который условно можно назвать итерациями Гаусса-Зейделя. Одна глобальная (5 + 1)-ая итерация заключалась в том, что в расчетной
_МЕЖДУНАРОДНЫЙ НАУЧНЫЙ ЖУРНАЛ «СИМВОЛ НАУКИ» №10/2015 ISSN 2410-700Х_
области Q последовательно перебирались все ячейки. При этом в каждой ячейке СЛАУ (14) решалось
ортогональным методом (Гивенса или Хаусхолдера), в правой части уравнений (11), (12) в качестве U , p берутся значения, известные на момент построения решения на (5 + 1)-ой итерации в данной ячейке.
Замечание. Пусть область Q разрезана на отдельные куски (частичные области), содержащие произвольное количество ячеек, линиями, проходящими по границам ячеек. Очевидно, что при итерационном решении дискретной задачи итерации в каждой частичной области можно вести параллельно во времени. При этом после одной (или нескольких) итераций перед следующей в правые части условий согласования (11), (12) для ячеек, которые попали на границу частичной области, нужно подставить
соответствующие, сосчитанные значения U , p из соседних к ним в Q ячеек, находящихся на границе
соседней частичной области. Здесь нет необходимости для обмена между решениями задачи в соседних частичных областях применять приближенные численные процедуры типа экстраполяции, которые приводят к потере достигнутой точности решения на текущей итерации и замедляют итерационный процесс решения глобальной задачи. Такой простой способ распараллеливания, который по сути заложен в формулах метода, был реализован в [10] на персональной ЭВМ с четырьмя процессорами и дал ускорение решения задачи примерно в 2.5 раза. С увеличением числа процессоров при распараллеливании решения задачи коэффициент ускорения должен возрастать.
2. Вариант алгоритма Крылова с редукцией базиса подпространства
Приведем кратко сведения о предложенном и более подробно описанном в [10] варианте алгоритма метода Крылова ускорения итерационного процесса решения СЛАУ. Пусть СЛАУ имеет вид
Х = Т'X + f, (15)
где вектор X — искомое решение, Т — в общем случае прямоугольная матрица с числом строк 6l' лыпим или равным числу неизвестных — компонент вектора X . / — вектор-столбец. Пусть матрица Т — матрица полного ранга, и пусть сходится итерационный процесс
Xn+l=TX"+f, п = ОД,..., (16)
в котором X" — приближение к решению на итерации с номером П . По определению, Г " = ТХ" + / — X" = Х"+] — X" — невязка уравнений (15), а /Г" = X — X" — погрешность
решения на fl -ой итерации. Нетрудно из выписанных формул вывести соотношение Т = ТГ ". Численное значение невязки на любой итерации в отличие от погрешности легко находится. Пусть от некоторого — начального приближения к решению сделаны k +1 итерации, то есть, вычислены
у-1 у-2 ук+1 -fO -1 ~к
Л ,Л ... Л и Г , Г ,...,/ , соответственно, о алгоритмах ускорения итерационных процессов по Крылову значение
до вычисления по формуле (16) уточняется, полагая
у*к+1 _ ук+1 ук+1
X — X + Y . Поправка вида
к
f*+1 =£«,/■' (17)
i—1
с неопределенными коэффициентами OCi ищется в подпространстве [4, 7]
К* (г1, Т) = span{r1, Т г1,..., Т*"1
Г }, где span{v ,..., V } — линейная оболочка векторов V , ...,
— к
V . Коэффициенты ОСх, ..., ОСк находятся из условия минимизации функционала невязки F(ai,...,ak) =
=\\Х к+1-ТХ k+l-f\\22
, которая возникает при подстановке
в (15). Здесь || U ||2 - евклидова
норма вектора U размерности N :
N
й\\1 =
2 i=1
2Х' u=(ul,...,uNу
Можно показать [6] , что для выполнения этого условия искомые коэффициенты ( разложения поправки в базисе подпространства должны удовлетворять следующей СЛАУ:
(У -r0)^ +... + {rk -fk~^ak =~fk. (18)
Построенное таким образом приближенное решение можно использовать как начальное
приближение для следующих итераций (16).
В случае сходящегося итерационного процесса Нщ " ||2= 0 . В малой окрестности решения
И—»со
элементы матрицы системы (18) — разности близких между собой малых величин г' — г' \ Z = 1,2,...Д. Из-за ошибок округления на ЭВМ их вычисление приводят к относительно большой погрешности. В результате по мере достижения всё большей точности решения СЛАУ поправки (17) вычисляются всё менее точно, рассматриваемый метод ускорения становится всё менее "устойчивым" и менее эффективным. Другая неприятность, которая снижает эффективность поправки в области малых невязок, заключается в том, что в условиях ограниченной разрядности представления чисел на ЭВМ система (18) нередко может быть плохо обусловленной или вырожденной. У нее с точностью до небольшого числа, близкого к машинному нулю, могут быть нулевые или линейно зависимые между собой столбцы.
Сделаем предположение (CR): допустим, что матрица системы (18) является матрицей полного ранга (complete rank). Для достижения большей устойчивости процесса вычисления коэффициентов
поправки Yl:+] в области малых невязок здесь применяются несколько приемов. Первый прием — нормировка столбцов матрицы системы (18). Это позволяет уменьшить в процессе решения СЛАУ (18) количество арифметических действий с числами, близкими к машинному нулю. Нормировка уравнений в
(18) производилась здесь при помощи замены неизвестных = OCj ||г' — Г7 ' 112. / = 1,...,АГ . В результате СЛАУ (18) принимает вид
Вр = (В1)р1+... + (Вк)рк=-г\ (19)
где Bi = (г1 —Г1 1)! | \г1 —Г1 1 112, i = 1,... - столбцы матрицы В .
Вторым приемом является применение ортогонального метода для решения переопределенной СЛАУ
(19). Здесь использовались методы либо Хаусхолдера, либо Гивенса с выбором главного элемента в i -ом столбце Вt. Z = 1,2,... Д. То есть, для прямоугольной матрицы В системы (19) строилось QR-разложение с ортогональной матрицей Q и прямоугольной матрицей R , у которой на главной диагонали ненулевые элементы, а под главной диагональю все элементы нули. Так что в первых к строках матрицы
R квадратная верхне-треугольная матрица. И решение переопределенной системы ВР=Гк (19) свелось к решению системы
Rp = -QTrk. (20)
Решение СЛАУ, состоящей из первых k уравнений системы (20) принимается в качестве псевдорешения системы (19). Известно, что ортогональный метод решение заданной системы сводит к решению системы, обусловленность которой не хуже обусловленности исходной. Описанные первый и второй приемы содержатся в упомянутом методе GMRES.
При численном решении прикладных задач в начале итерационного процесса, когда невязки Г1 (i = 0,1,...) не маленькие, с большой вероятностью матрица В имеет полный ранг — предположение
CR выполняется. Однако по мере сходимости итерационного процесса при небольших Г с ненулевой
вероятностью предположение CR может не выполниться. При этом попытка построить QR -разложение
не дает в первых к строках матрицы R треугольную матрицу размера к . Если этот факт имеет место с точностью близкой к машинному нулю, т.е. по крайней мере один из элементов на диагонали верхней треугольной матрицы R близок к машинному нулю, то поправка (17) будет найдена с непригодной точностью (или попытка решить систему (20)) может дать в результате большое число, не представимое в арифметике с конкретной используемой разрядностью представления чисел в ЭВМ). Это может привести (или приведет) к автоматическому останову машины (авосту). Если даже система (20) будет решена, то это
может быть сделано с плохой точностью.При этом поправки Yl:+] к Xk+l вычисляются неустойчиво. В вычислительном эксперименте в области малых невязок решения можно наблюдать осцилляции в графике погрешности решения.
Третий существенный прием, примененный здесь, как и первые два направлен на увеличение устойчивости вычисления поправок (устойчивого построения подпространства Крылова). Суть его
заключается в редукции имеющегося базиса путем ограничения числа невязок Т1, i = 0,1,..S < к , используемых в редуцированном базисе, при вычислении поправок к приближенным значениям решения в области плохой обусловленности системы (18). Здесь применялся достаточно простой критерий (S -критерий) для выбора числа невязок, пригодных с точки зрения устойчивого построения подпространства Крылова. В [7] был реализован некоторый его вариант.
Пусть в процессе QR -разложения на диагонали матрицы R в (s + 1) -ой строке появился элемент
по модулю меньше некоторого числа 8, близкого к машинному нулю, и все элементы под диагональю равны нулю. Тогда для построения поправки используются только первые S +1 невязка Т1 (J = 0,1,.. .,5"). В системе (19) оставляются первые S столбцов, а в правой части соответственно
к s
редуцированной определяющей системы вместо вектора —Г берется вектор —Г .
Как показали результаты большого количества итерационного решения тестовых примеров, применение S-критерия позволяет избежать некоторого количества авостов при приближении величин невязок к машинному нулю и получить более точные решения СЛАУ, чем без применения критерия. Очевидно, что программа с автоматическим применением S -критерия позволяет применять метод Крылова также в случае, когда матрица B системы (19) не является матрицей полного ранга. 3. Ускорение сходимости итераций применением
многосеточного алгоритма Классическим способом реализации многосеточных алгоритмов являются V- и W-циклы [9]. В V-цикле расчет начинается на самой мелкой пространственной сетке (сетке с самыми мелкими ячейками)
1 - К ,
многосеточного комплекса, имеющей шаг = , где К - шаг самой грубой сетки, m - количество
сеток, используемых в V-цикле, m >1. Затем на нисходящей ветви V-цикла осуществляются расчеты на
7 - К
последовательности все более грубых сеток, при этом переход с мелкой сетки с шагом К _
j j'
__К
т
сужения (restriction).
] — \ . — 1 на более грубую сетку с шагом /?/+| — —т-]-\ осУЩествляется с помощью оператора
_МЕЖДУНАРОДНЫЙ НАУЧНЫЙ ЖУРНАЛ «СИМВОЛ НАУКИ» №10/2015 ISSN 2410-700Х_
После расчета на самой грубой сетке, начинается движение операциями продолжения по восходящей ветви V-цикла — расчеты на последовательности все более мелких сеток. Переход с более грубой сетки на более мелкую осуществляется с помощью операторов продолжения (prolongation).
В методе КНН, как и в других методах, количество итераций, необходимых для достижения заданной точности приближения к предельному для конкретного итерационного процесса решению зависит от начального приближения: чем ближе начальное приближение к предельному решению, тем меньше итераций требуется для достижения заданной точности приближения к нему. Как способ получения хорошего начального приближения для итераций на самой мелкой сетке из числа используемых в многосеточном комплексе, сначала на самой грубой сетке задается нулевое начальное приближение для
коэффициентов разложения b j т в каждой ячейке сетки; затем решается СЛАУ на этой грубой сетке (что
требует небольшого машинного времени). Далее полученное решение используется в качестве начального приближения на следующей более мелкой сетке. Продолжая этот процесс на последовательности все более мелких сеток, получаем хорошее начальное приближение решения для итераций на самой мелкой сетке.
В случае метода КНН операцию продолжения можно осуществить без потери точности построенного на грубой сетке решения, избегая интерполяции и любых приближенных численных операций. Для этого достаточно в каждой ячейке более мелкой сетки переписать в ее локальных координатах полиномиальное решение, полученное в ячейке грубой сетки, содержащей эту ячейку более мелкой сетки.
Проиллюстрируем алгоритм операции продолжения на примере составляющей скорости Щ (y,
У2,Ь1,..., Ь]2) . Пусть hx = /?. где H — полушаг грубой сетки, и пусть /?2 = /?, /2 — полушаг мелкой
сетки, на которой нужно найти разложение функции u по базису.
Шаг 1. Пусть X1, X2 — глобальные координаты центра ячейки грубой сетки. Сделаем следующие
подстановки в полиномиальное выражение для Щ :
y =(xl-Xl )/ h, l = 1,2. (21)
В результате получим полином
Ul(x1,x2,bl,...,bl2) = ul ——-Д,...,^
12
(22)
V "1 "1 /
Шаг 2. Пусть (Х], X 2 ) — глобальные координаты центра любой из четырех ячеек мелкой сетки,
содержащихся в ячейке грубой сетки. Сделаем в (22) замену х1 = Х/ + у/ ■ [%. / = 1,2. В результате
получим многочлен второй степени = Р(ух,у2,Ьх,.. от переменных у ,у2 -коэффициентами
,..., 2 . После приведения подобных оказывается, что координаты Х], X2 и Х^, X2 входят в Ь,
(/ = 1,...,12) только в виде комбинаций 8х1 — (Х1 — Х;) / . Согласно (21), величина
~8х1 = (X! — X!) / является локальной координатой в ячейке грубой сетки центра ячейки мелкой сетки.
Приведем выражения для коэффициентов Ь (7 = 1,...,12) представления решения в ячейке мелкой сетки с полушагом Н2 в терминах коэффициентов Ьх,..., Ь]2 представления решения в ячейке с полушагом Ь = 2к2 и величин 71 = Ъ2 / ^, 72 = (( :
ЪХ=ЪХ- Ъ28хх + Ъ48х\ - (Ъъ + 2 Ь58хх)8х2 + Ь68х2; Ь2 = <зх(Ь2 - 2 Ъ48хх + 2Ь58х2); 4 = сг, [Ъъ + 2(Ь58х] - Ь68х2 )]; Ъ4 = сг2Ь4; Ь5 = сг2Ь5; Ь6 = сг2Ь6;
b1=b1- b%8x] + b9Sx\ + Sx2(b2 - 2b4Sxl + b5Sx2); bH = cr,(/)8 - 2b9Sxx + 2b4Sx2);
2/'
b9-cr2b9; bw-bw
■ Ъх х5хх - Ъх2дх2; 1 = 1; Ьп= ■ Нетрудно видеть, что можно выписать аналогичные формулы перехода от грубой сетки к мелкой в случае любой целой величины отношения сторон их ячеек. Обратно: пусть требуется начать итерации на грубой сетке с использованием полученного решения на более мелкой сетке после 5-ой итерации. Заметим, что в коэффициенты и правые части уравнений (14) для вычисления решения на 5 + 1-ой итерации входят значения решения с предыдущей 5-ой итерации. Для начала счета на грубой сетке в коэффициенты и правые части уравнений (14) надо подставить значения решения, полученного на мелкой сетке, определяя по формулам замены координат, какие точки в ячейках мелкой сетки соответствуют точкам коллокаций и согласования решения, участвующим в записи уравнений (14) в ячейке грубой сетки.
Таким образом, при переходе с сетки на сетку в операциях продолжения и сужения на многосеточном комплексе в методе КНН достаточно выполнить точную операцию — замену координат, которая в отличие от таких приближенных численных операций, как интерполяция и другие, не вносит дополнительной (своей) погрешности в построенное на промежуточной итерации решение и тем самым не сказывается отрицательно на скорости сходимости решения. В этом одно из положительных проявлений наличия кусочно--аналитического решения дискретной задачи и специфика при использовании многосеточных вариантов метода КНН.
4. Результаты численных экспериментов 4.1. Тестирование
Рассмотрим следующее точное решение уравнений Навье-Стокса (1), (2) в квадратной области (5):
щ = соэ(2жх! )эт(27х2), и2 = -эт^ж^ )соэ(27Х2),
p
соэ
7ХЛ
+ соэ
7Х
2 X sin(^X /2)
(23)
Ж
Заметим, что решение (23) удовлетворяет уравнению неразрывности (2). Выражения для правых части , уравнений (1) легко получаются подстановкой решения (23) в (1), поэтому здесь они не
приводятся. Тестовое решение (23) не имеет особенностей, поэтому можно наблюдать порядок сходимости численного решения уже на грубых сетках. Область (5) покрывалась равномерной сеткой квадратных ячеек. Припишем индексы I, ] ячейки к ее геометрическому центру, где I, ] меняются, соответственно,
вдоль осей Х, Х2 .
Поскольку описываемый метод использует линеаризацию уравнений Навье-Стокса, необходимы
итерации по нелинейности. Использовалось нулевое начальное приближение для величин Щ и р .
Для того, чтобы определить численные значения погрешности решения на конкретной равномерной
сетке с полушагом И вычислялись следующие среднеквадратичные величины:
-|0.5
1 ММ 2
-ИЖ,,] - *
Su(h) =
2M1
u
)2
i=1 j=1 V=1
, Sp(h) =
J MM
11 i=1 J=1
Pi'j)
0.5
Мех ех
— количество ячеек вдоль каждого координатного направления, Щ . и р J — вектор скорости и давление, вычисленные из точного решения (23). Величины Щ и р обозначают численное решение,
полученное по методу КНН, описанному выше. Будем вычислять порядки сходимости Уи и V из численных решений для вектора скорости и и для давления р по следующим формулам, известным в численном анализе:
_МЕЖДУНАРОДНЫЙ НАУЧНЫЙ ЖУРНАЛ «СИМВОЛ НАУКИ» №10/2015 ISSN 2410-700Х_
v = \og[Su(hm_l)]-log[Su(hm)] ^ = \og[Sp(hm_,)]-\og[Sp(hm)]
" 10g(hm_1)-10g(hm) ' * lOg(Ä^l)- 10g(hm) '
где hm, 171= 2,3,... — некоторые значения шага И. такие, что Ф Ит. Пусть Ь '^ . /.
S = 0 1 ... — значение коэффициента Ъ. , в (9) на S -ой итерации. Использовалось следующее условие для окончания итераций по нелинейности: Sbs+1 < £ , где
Sbs+l = max(max Ц+^^Ц ^ ), (25)
i j \1</ <12 /
и £ < h2 — малая положительная величина. В дальнейшем будем называть величину (25) псевдопогрешностью приближенного решения.
Изучалось влияние числа уравнений коллокации и его соотношения с числом уравнений согласования на точность расчетов. Сравнительно хорошая точность решения получалась, когда число уравнений коллокаций примерно равнялось или превосходило число уравнений согласования не более, чем в два раза. Надо также учитывать, что величина погрешности округлений в решении задачи связана с числом уравнений в СЛАУ дискретной задачи. И если число уравнений в СЛАУ дискретной задачи превосходит некоторую величину, то погрешность округлений в решении ухудшает точность и сходимость решения.
1000 2000 3000 4000 5000 6000 0 1000 2000 3000 4000 5000 6000 (а) (б)
Рисунок 2 - Влияние значения к в (18) на скорость сходимости величин log (fib" (рис. (а)) и
log fiu (рис. (б)), где n — число итераций.
Рис. 2 иллюстрирует влияние величины к в (18) на скорость сходимости метода КНН. Эти расчеты были выполнены на сетке из 20 X 20 ячеек. Число Рейнольдса К_е=1000, X = 0.5 в (5); 77 =3.5,
72 =1.0 в (11). Критерием останова счета было выполнение неравенства ЗЪ" <10 9 . Расчету без
применения алгоритма Крылова соответствует случай к = 0 . В этой серии расчетов в переопределенную СЛАУ (14) было включено уравнение (13). В процессе итераций по методу КНН абсолютная величина
интеграла (4) падала со значения порядка 10 3 до величины порядка 10 12 —10 13, то есть до величины порядка машинных ошибок округления при расчетах по Фортран-программе с двойной точностью. Это может служить одним из критериев правильности программной реализации представленного выше метода КНН. Видно, что с увеличением числа невязок к в методе Крылова скорость сходимости численного решения по методу КНН растет.
Количество итераций Ыц, необходимых для достижения неравенства ЗЪ" <10—9 , составляло 29115, 1852, 1233 и 1072, соответственно, при к = 0, 2, 10 и 20. Таким образом, применение алгоритма Крылова при
_МЕЖДУНАРОДНЫЙ НАУЧНЫЙ ЖУРНАЛ «СИМВОЛ НАУКИ» №10/2015 ISSN 2410-700Х_
Re = 1000 с k = 20 позволило уменьшить количество итераций, требуемых для сходимости приближенного решения, в 27.16 раз
Таблица 2
Погрешности Su,Sp и их порядки сходимости Vu, Vp на последовательности сеток, Re = 1000,
X = 0.25
Su Sp Vu VP
20 0.342 -10-3 0.262 -10-3
40 0.982 -10-4 0.144 -10-3 1.80 0.86
60 0.412 -10-4 0.779 -10-4 2.14 1.51
80 0.161-10-4 0.486 -10-4 3.26 1.65
96 0.108 -10-4 0.347 -10-4 2.19 1.84
по сравнению со случаем к = 0. Аналогичная динамика уменьшения N при использовании алгоритма Крылова была получена и при Re = 100.
На рис. 2 видно, что, начиная примерно с к = 10, результаты ускорения итераций по Крылову незначительно различаются между собой по количеству затраченных итераций N . Уравнения Навье-
Стокса при умеренных и больших числах Рейнольдса существенно нелинейны. Значения решения уравнений "сильно" зависят от значений ее коэффициентов. Решения нелинейных уравнений только в достаточно малой окрестности решения (назовем его опорным), на котором линеаризованы уравнения, мало отличаются от решения линеаризованных уравнений. Поправленное решение с увеличением к все больше и больше удаляется от опорного решения. И вне некоторой достаточно малой окрестности опорного решения поправки становятся менее эффективными. Кроме того, с увеличением числа векторов невязок к в системе (18) увеличиваются необходимые для ее решения число арифметических действий и машинное время, растет величина накапливаемой ошибки округлений, и уменьшается точность поправок. Поэтому во многих случаях применения вышеописанного метода КНН в сочетании с алгоритмом Крылова для ускорения итераций при решении уравнений Навье-Стокса с умеренными числами Рейнольдса целесообразно брать количество невязок к из интервала 2 < к < 10.
Из табл. 2 видно, что при Re = 1000 порядок сходимости V лежит в интервале 1.80 < V < 3.26 , и
при М > 60 он выше, чем второй порядок. При этом V < V для всех рассмотренных сеток.
Были проведены расчеты с применением восходящей ветви У-цикла с целью выяснить, как влияет применение только многосеточного алгоритма на ускорение сходимости метода КНН. Были также проведены расчеты, в которых движение по восходящей ветви многосеточного У -цикла применялось совместно с алгоритмом ускорения, основанным на подпространствах Крылова. Результаты этих расчетов представлены в табл. 3. В ней Кт§г — количество последовательно используемых сеток в многосеточном
комплексе. Если =1, то это означает, что в расчете используется только одна сетка, и это самая
мелкая сетка с числом ячеек 80 X 80. Далее, — суммарное количество итераций, выполненных на всех сетках комплекса. Фактор ускорения итерационного процесса AF в результате применения того или иного способа его ускорения вычисляется как отношение времени счета при =1, к = 0 ко времени
счета при применении последовательности из нескольких сеток (>1, к = 0) или же последовательности сеток в сочетании с применением алгоритма Крылова на каждой сетке
(Кт§г >1, к >1).
Критерием останова счета было выполнение неравенства Sbn <10"9 . Значения погрешностей Su и Sp даны на момент завершения расчета. На рис. 3 показан фактор ускорения AF как функция числа ячеек
M flne х M ßne, где M ßne - число ячеек самой мелкой сетки вдоль каждого координатного
направления, которая используется в конкретной последовательности сеток. Для построения графика AF были взяты самые высокие факторы ускорения, полученные на различном количестве сеток,
использовавшихся в многосеточном алгоритме. Видно, что с ростом M ßne в приведенном диапазоне имеет место нелинейный рост фактора ускорения AF.
Таблица 3
Влияние комбинации методов Крылова и Федоренко на последовательности сеток размеров 5-2'"х5-2'", т = 0,...,4 на скорость сходимости метода КНН при Re = 1000.
Метод N Tсчета ' сек. AF Su Sp
Kmgr =1, k = 0 3827643 697397 1.0 0.2801 -10 4 0.5428 -10~3
Kmgr = 5, к = 0 219375 40904 17.05 0.5966 -10~4 0.6314 -10~4
Kmgr =1, к = 9 39871 7454 93.56 0.9310 -10~4 0.1795 -10~3
K =1, к = 10 mgr ? 40278 7544 92.44 0.9310 -10~4 0.1796 -10~3
K =5, к = 5 mgr ? 16685 2856 244.19 0.5976 -10~4 0.9101-10 4
Kmgr =5, к = 7 14171 2398 290.82 0.5966 -10 4 0.1135 -10~3
Kmgr =5, к = 8 13730 2318 300.86 0.5976 -10~4 0.1086 -10~3
K =5, к = 9 mgr ? 12851 2157 323.32 0.5976 -104 0.1077 -103
Kmgr =5, к = 10 15552 2676 260.64 0.5976 -10^ 0.9164 -10~4
Рис. 4 иллюстрирует влияние применения последовательности сеток и алгоритма Крылова на скорость сходимости метода КНН (последовательность сеток и число Рейнольдса те же, что и в случае табл. 2). Эти
300 250 200 J50 100 50
£F
%1пс
20
40
60
Ж)
Рисунок 3 - Фактор ускорения AF как функция числа ячеек M jine X M самой мелкой сетки, использованной в методе Федоренко.
Рисунок 4 -Результаты по ускорению решения тестовой задачи методом КНН применением методов Крылова и Федоренко при Re = 1000. На
оси абсцисс отложена величина + 1) ,
где Т — время решения задачи в минутах.
_МЕЖДУНАРОДНЫЙ НАУЧНЫЙ ЖУРНАЛ «СИМВОЛ НАУКИ» №10/2015 ISSN 2410-700Х_
Таблица 4
Результаты по ускорению решения методом КНН эталонной задачи путем применения комбинации методов Крылова и Федоренко на последовательности сеток размеров 5 • 2m х 5 • 2m, m = 0,1 ,2,3 при
Re = 100
Метод N , сек. AF
Kmgr =1 k = 0 3377640 139598 1.0
Kmgr =4, k = 0 166501 6347 21.99
Kmgr =1, k = 10 137880 5860 23.82
Kmgr = 4, k = 8 7593 277 503.96
Kmgr =4, k = 9 6121 220 634.54
Kmgr =4, k = 10 6524 238 586.55
численные эксперименты показывают, что с увеличением размера сеток (до некоторого предела) и требуемой точности решения итерационный процесс, основанный на указанной комбинации многосеточного алгоритма с методом Крылова, становится все более эффективным в сравнении с процессом на самой мелкой сетке многосеточного комплекса. Это четко отражено в табл. 2 и на рис. 3. 4.2. Использование высокоточных решений эталонной задачи
для верификации численных экспериментов Для контроля и верификации результатов численных экспериментов здесь кроме их сравнения с точными решениями тестовых задач проведено сравнение характерных величин, полученных в численных экспериментах по решению эталонной задачи о течении вязкой жидкости в каверне. Результаты все возрастающей точности ее решения публикуются различными исследователями на протяжении последних тридцати лет в подробных таблицах. Наиболее точные решения этой задачи (на конкретный текущий момент времени) [7, 21, 27] можно использовать для выяснения точности новых и модифицированных существующих численных методов решения уравнений Навье-Стокса. В обсуждаемой задаче
1.0
х2 0.8
0.6
0.4
0.2
(а) (б)
Рисунок 5 - Картина линий тока (а) и профиль составляющей скорости VI вдоль линии х\ = 0.5 (б) в решении эталонной задачи при Re = 1000 (значки А — результаты работы [27]).
_МЕЖДУНАРОДНЫЙ НАУЧНЫЙ ЖУРНАЛ «СИМВОЛ НАУКИ» №10/2015 ISSN 2410-700Х_
расчетная область каверна - квадрат (5) со стороной X = 1, начало координат находится в ее левом нижнем углу. Верхняя крышка каверны движется в безразмерных величинах с единичной скоростью в
положительном направлении оси Ox^. Остальные стороны каверны (5) покоятся. На всех сторонах заданы
условия прилипания: V =1, V =0 при Х2 = X и V =0 , m = 1,2, на остальных сторонах.
Течение в каверне с движущейся крышкой имеет сингулярности в верхних углах области. Их влияние на точность численного решения усиливается с увеличением числа Рейнольдса. Поэтому при больших числа Рейнольдса для получения более точного решения необходимо применять адаптивные сетки с более мелкими ячейками в окрестности сингулярностей. Здесь использовали только равномерные сетки, размер которых не превосходил 320 X 320 ячеек.
Таблица 4 иллюстрирует эффективность совместного применения алгоритмов Крылова и Федоренко.
Видно, что ускорение сходимости весьма значительное: при Kmgr =4, k = 9 соответствующий фактор
ускорения AF = 634.54. Из этой таблицы также следует, что в рассматриваемой задаче наилучшим для ускорения сходимости итерационного процесса при решении уравнений Навье-Стокса с
умеренными числами Рейнольдса является число невязок k = 9 в описанном выше варианте алгоритма Крылова. Заметим, что это же значение k оказалось наилучшим и при тестировании совместного применения алгоритмов Крылова и Федоренко на решении (23), см. табл. 2.
На рис. 5 представлены некоторые результаты численных расчетов течения вязкой несжимаемой жидкости в квадратной каверне для числа Рейнольдса Re = 1000. Расчеты выполнены на равномерной сетке 320 X 320 ячеек. Стрелки на рис. 4, (а) указывают локальные направления движения частиц жидкости. Сплошная линия на рис. 4, (б) — результат расчета по вышеописанному методу КНН. 4.3. Влияние свойств метода КНН на эффективность алгоритмов
ускорения итерационных процессов Достигнутое существенное ускорение сходимости итераций при комбинированном применении методов Крылова и Федоренко в методе КНН в сравнении с аналогичным способом ускорения итераций в методе конечных объемов [14] обусловлено несколькими причинами. Во-первых, в методе КНН при переходе с сетки на сетку на многосеточном комплексе в качестве начального приближения на текущей сетке берется решение на последней итерации на предыдущей сетке, в котором делается, как показано выше, только точное преобразование — замена координат (см. п. 3). То есть, не используются приближенные численные операции типа интерполяции (экстраполяции, сглаживания), которые вносят в решение на промежуточном этапе дополнительную погрешность, ухудшают достигнутую в нем промежуточную точность и тем самым замедляют сходимость итераций.
Во-вторых, при применении ортогональных методов решения переопределенных СЛАУ не ухудшается обусловленность дискретной задачи в отличие от применения для этого собственно МНК. Применение МНК в перечисленных в обзоре во введении случаях в отличие от применения ортогональных методов в силу ухудшения обусловленности СЛАУ замедляет сходимость итераций и снижает эффективность от применения методов ускорения итераций.
В третьих, в методах КНН и КНК минимизируется функционал невязки. Тем самым на каждой промежуточной итерации независимо от применения подпространств Крылова и многосеточного комплекса подавляются невязки уравнений дискретной задачи и связанные с ними погрешности численного решения. Это также независимо от других указанных факторов вносит свой вклад в ускорение сходимости итераций. При этом в методе КНН прежде всего в задачах с большими СЛАУ влияние этого фактора на итерационный процесс их решения по указанным выше причинам сильнее, чем в методе КНК. В силу перечисленных причин в первом случае больше возможностей решить задачу на ЭВМ, чем во втором.
Для иллюстрации роли минимизации невязки уравнений дискретной задачи в методах КНН и МНК на рис. 6 приведен типичный результат решения методом наименьших квадратов задачи ПАФ с некоторой погрешностью во входных данных, имеющей случайный характер. В различных экспериментах брались различные гладкие тестовые функции f(x) (конкретно на рис. 6 f (x) = exp( x)). Их
/
Рисунок 6 - Сплошная линия - решение задачи ПАФ, штриховая линия - тестовая функция, • - входные данные задачи, полученные случайными возмущениями значений тестовой функции в узлах равномерной сетки.
дискретные значения в узлах равномерной сетки по координате X в интервале (0, 1), возмущенные датчиком случайных чисел, брались в качестве входных данных задачи ПАФ. В этих экспериментах в качестве невязки бралась разница значений координаты по оси ординат точек входных данных и соответствующих им в узлах сетки точек аппроксиманта, построенного из требования минимизации функционала невязки (сумма квадратов невязок во всех узлах сетки). Нетрудно видеть, что в узлах сетки точки на аппроксиманте расположены ближе к тестовой функции, чем соответствующие им точки входных данных задачи. Более того, ни один из полиномов первой, второй и других степеней в качестве лагранжевого интерполянта не может приблизить точнее тестовую функцию на всем интервале (0, 1) по указанным входным данным, чем приведенное на рисунке решение задачи ПАФ. Имеет место аналогия между решением задачи ПАФ методом наименьших квадратов и решением путем минимизации функционала невязки уравнений дискретной задачи в методах КНК и КНН.
5. Заключение
Представлено несколько новых вариантов алгоритмов ускорения сходимости итераций при численном решении СЛАУ, возникающих в результате дискретизации PDE методом КНН. Один из них — вариант метода Крылова с редукцией базиса подпространства Крылова в области малых невязок итерационнго решения СЛАУ. Применение редукции базиса позволило в этой области более точно строить поправки по Крылову решение СЛАУ, чем в случае неприменения редукции. Другой алгоритм ускорения итерационного решения дискретной задачи в методе КНН комбинирует операцию продолжения метода Федоренко с новым вариантом метода Крылова. Проведены вычислительные эксперименты с комбинированным применением указанных двух методов, а также с их раздельным применением при численном решении уравнений Навье-Стокса. Наибольший эффект ускорения сходимости рассмотренных итерационных процессов дал комбинированный алгоритм. При этом достигнут выигрыш во времени решения дискретной задачи в сотни раз по сравнению со случаем ее решения лишь на самой мелкой сетке использованного многосеточного комплекса без применения методов Крылова и Федоренко. Отмечено, что наличие кусочно - аналитического решения дискретной задачи в методе КНН, минимизация функционала невязки уравнений дискретной задачи и применение ортогональных методов решения СЛАУ положительно влияют на эффективность применения такого способа ускорения итерационного решения задачи. Рассмотренные алгоритмы легко могут быть реализованы при решении трехмерных задач и при применении других численных итерационных методов решения PDE.
Список использованной литературы
1. Edwards W.S., Tuckerman L.S., Friesner R.A., Sorensen D.C. Krylov methods for the incompressible Navier--Stokes equations // J. Comp. Phys. 1994. V. 110. P. 82-102.
2. Knoll D.A., Keyes D.E. Jacobian-free Newton--Krylov methods: a survey of approaches and applications // J. Comp. Phys. 2004. V. 193. 2. P. 357-397.
3. Griffith B.E. An accurate and efficient method for the incompressible Navier-Stokes equations using the projection method as a preconditioner // J. Comp. Phys. 2009. V. 228. No. 20. P. 7565-7595.
4. Saad Y. Numerical methods for large eigenvalue problems. Manchester: Manchester University Press, 1991.
_МЕЖДУНАРОДНЫЙ НАУЧНЫЙ ЖУРНАЛ «СИМВОЛ НАУКИ» №10/2015 ISSN 2410-700Х_
5. Крылов А.Н. О численном решении уравнения, которым в технических вопросах определяются частоты малых колебаний материальных систем // Изв. АН СССР, Отд. матем. и естеств. наук. 1931. № 4. С. 491539.
6. Слепцов А.Г. Об ускорении сходимости линейных итераций II // Моделирование в механике. Новосибирск, 1989. Т. 3(20). № 5. C. 118-125.
7. Исаев В.И., Шапеев В.П. Варианты метода коллокаций и наименьших квадратов повышенной точности для численного решения уравнений Навье-Стокса // Ж. вычисл. матем. и матем. физ. 2010. Т. 50. № 10. С. 1758-1770.
8. Федоренко Р.П. О скорости сходимости одного итерационного процесса. // Ж. вычисл. матем. и матем. физ. 1964. Т. 4. № 3. С. 559-564.
9. Wesseling P. An Introduction to Multigrid Methods. Chichester: Wiley & Sons, 1992.
10. Шапеев В.П., Ворожцов Е.В., Исаев В.И., Идимешев С.В. Метод коллокаций и наименьших невязок для трехмерных уравнений Навье--Стокса // Вычислит. методы и программирование. 2013. Т. 14. С. 306-322.
11. Piquet, J., Vasseur X. A nonstandard multigrid method with flexible multiple semicoarsening for the numerical solution of the pressure equation in a Navier-Stokes solver // Num. Algorithms. 2000. V. 24. No. 4. P. 333-355
12. Jothiprasad G., Mavriplis D.J., Caughey D.A. Higher-order time integration schemes for the unsteady Navier-Stokes equations on unstructured meshes // J. Comp. Phys. 2003. V. 191. No. 2. P. 542-566.
13. Sotiropoulos F. A numerical method for solving the 3D unsteady incompressible Navier-Stokes equations in curvilinear domains with complex immersed boundaries // J. Comp. Phys. 2007. V. 225. No. 2. P. 1782-1809.
14. Lucas P., van Zuijlen A.H., Bijl H. Fast unsteady flow computations with a Jacobian-free Newton-Krylov algorithm // J. Comp. Phys. 2010. V. 229. No. 2. P. 9201-9215.
15. Nasr-Azadani M.M., Meiburg E. TURBINS: An immersed boundary, Navier-Stokes code for the simulation of gravity and turbidity currents interacting with complex topographies // Comp. & Fluids. 2011. V. 45. No. 1. P. 1428.
16. Benzi M., Wang Z. Analysis of augmented lagrangian-based preconditioners for the steady incompressible navier-stokes equations // SIAM J. Sci. Comput. 2011. V. 33. No. 5. P. 2761-2784.
17. Fairag F.A., Wathen A.J. A block preconditioning technique for the streamfunction-vorticity formulation of the Navier-Stokes equations // Num. Methods Partial Differential Equations. 2012. V. 28. No. 3. P. 888-898.
18. Wang M., Chen L. Multigrid methods for the stokes equations using distributive gauss-seidel relaxations based on the least squares commutator // J. Sci. Comput. 2013. Vol. 56. No. 2. P. 409-431.
19. Nickaeen M., Ouazzi A., Turek S. Newton multigrid least-squares FEM for the V-V-P formulation of the Navier-Stokes equations // J. Comp. Phys. 2014. V. 256. P. 416-427.
20. Беляев В.В., Шапеев В.П. Метод коллокаций и наименьших квадратов на адаптивных сетках в области с криволинейной границей // Вычисл. технологии. 2000. Т. 5. № 4. С. 12-21.
21. Shapeev A.V., Lin P. An asymptotic fitting finite element method with exponential mesh refinement for accurate computation of corner eddies in viscous flows // SIAM J. Sci. Comput. 2009. V. 31. No. 3. P. 1874-1900.
22. Shapeev V.P., Vorozhtsov E.V. CAS application to the construction of the collocations and least residuals method for the solution of the Burgers and Korteweg--de Vries-Burgers equations // Lecture Notes in Computer Science. Vol. 8660. Heidelberg: Springer, 2014. P. 432-446.
23. Плясунова А.В., Слепцов А.Г. Коллокационно-сеточный метод решения нелинейных параболических уравнений на подвижных сетках // Моделирование в механ. 1987. Т. 1(18). № 4. С. 116-137.
24. Jiang B.-N., Carey G.F. Adaptive refinement for least-squares finite elements with element-by-element conjugate gradient solution // Int. J. Numer. Meth. Engng. 1987. Vol. 24. P. 569-580.
25. Семин Л.Г., Слепцов А.Г., Шапеев В.П. Метод коллокаций-наименьших квадратов для уравнений Стокса // Вычисл. технологии. 1996. Т. 1. № 2. С. 90-98.
26. Исаев В.И., Шапеев В.П. Метод коллокаций и наименьших квадратов повышенной точности для решения уравнений Навье-Стокса // Докл. Акад. наук. 2012. Т. 442. № 4. С. 442-445.
_МЕЖДУНАРОДНЫЙ НАУЧНЫЙ ЖУРНАЛ «СИМВОЛ НАУКИ» №10/2015 ISSN 2410-700Х_
27. Botella O., Peyret R. Benchmark spectral results on the lid-driven cavity flow // Comput. & Fluids. 1998. V.
27. No. 4. P. 421-433.
28. Shapeev V.P., Vorozhtsov E.V. CAS application to the construction of the collocations and least residuals method for the solution of 3D Navier-Stokes equations // Lect. Notes Comp. Sci. Vol. 8136. Heidelberg: Springer, 2013. P. 381-392.
© Е.В. Ворожцов, В.П. Шапеев, 2015
УДК 004.42
Каменева Инга Олеговна
к.пед.н, доцент Каменев Владимир Викторович
к.физ.-мат., доцент ФМИТ МГУ г.Саранск, Российская Федерация [email protected]
ОБУЧЕНИЕ ОСНОВАМ ПРОГРАММИРОВАНИЯ. ОПТИМИЗАЦИЯ СКОРОСТИ РАБОТЫ
ПРОГРАММ
Аннотация
Рассматриваются приемы, позволяющие увеличить скорость выполнения программы, которым необходимо научить студентов при изучении основ программирования.
Ключевые слова Оптимизация, программа, скорость выполнения.
Программы, создаваемые студентами в процессе изучения дисциплин, связанных с основами программирования, не отличаются большими объемом и сложностью. Но, несмотря на это, преподаватель должен уже в начальной стадии овладения профессиональными навыками акцентировать внимание студентов на приемах, позволяющих создать эффективную программу.
Для пользователя основными факторами при оценке эффективности приложения являются скорость и размер. Вероятность одновременной оптимизации по обоим параметрам мала. Действительно, подход, позволяющий уменьшить размер, часто приводит к уменьшению скорости работы. И напротив -программа, оптимизированная по скорости, занимает больший объем, чем более «медленный» вариант. Преподаватель должен сказать и о том, что иногда изменения, сделанные для уменьшения размера или для увеличения скорости, приводят к ухудшению читаемости кода, такой код становится труднее поддерживать. Поэтому оптимизация - это не только овладение определенными методами, это в некотором роде искусство, т.к. программист должен уметь определять, где и когда следует применять эти методы.
Для максимальной оптимизации необходимо сначала определить так называемые «узкие места», т.е. фрагменты кода, которые являются основными потребителями ресурса. Если большую часть времени работы программы занимают вычисления, то для того, чтобы увеличить ее скорость студенты должны овладеть приемом, который называется понижением силы операций. Самые быстрые операции - это сложение и вычитание, медленнее выполняется умножение, еще медленнее выполняется деление. В
