close

Вход

Забыли?

вход по аккаунту

?

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

код для вставкиСкачать
Анализ эффективности итерационных
методов решения систем линейных
алгебраических уравнений,
реализованных в пакете OpenFOAM
И. К. Марчевский, В. В. Пузикова
МГТУ им. Н.Э. Баумана
iliamarchevsky@mail.ru, valeria.puzikova@gmail.com
Аннотация. Для выбора оптимального в смысле вычислительной эффективности
итерационного метода решения систем линейных алгебраических уравнений,
возникающих при дискретизации дифференциальных уравнений в частных
производных, помимо скорости сходимости следует учитывать такие характеристики
системы и метода, как число обусловленности, коэффициент сглаживания, показатель
«затратности». Последние две характеристики вычисляются по коэффициентам
усиления гармоник, которые позволяют судить о сглаживающих свойствах
итерационного метода и его «затратности», т.е. о том, насколько медленнее метод
подавляет низкочастотные компоненты ошибки по сравнению с высокочастотными.
Предложен способ определения коэффициентов усиления гармоник, основанный на
использовании дискретного преобразования Фурье. В качестве примера приведён
анализ эффективности метода BiCGStab c ILU и многосеточным предобуславливанием
при решении разностных аналогов уравнений Гельмгольца и Пуассона, возникающих
при моделировании течения вязкой несжимаемой жидкости в квадратной каверне.1
Ключевые слова. Разреженные линейные системы; предобуславливание;
сглаживатели; коэффициенты усиления гармоник; многосеточные методы.
1. Введение
Как правило, значительная доля всего объёма вычислительной работы при
численном моделировании технических систем, физических явлений и
технологических процессов приходится на решение систем линейных
алгебраических уравнений (СЛАУ), возникающих при дискретизации
соответствующих дифференциальных или интегро-дифференциальных
уравнений.
Существует несколько классов итерационных методов решения СЛАУ [1],
различающихся самим подходом к построению очередного итерационного
приближения. Это методы, основанные на расщеплении, а также методы
вариационного и проекционного типов. Необходимо отметить, что это
исторически сложившееся деление весьма условно, поскольку любой
итерационный метод можно записать в форме метода, основанного на
расщеплении.
Во многих программных пакетах, используемых при решении широкого
класса задач механики сплошной среды (ANSYS, NASTRAN, FLUENT, STARCD, COSMOSFloWorks и др.), имеется возможность выбора итерационного
метода, который будет использоваться при проведении расчетов. В некоторых
случаях методы решения линейных систем в программном пакете выбраны
«по умолчанию» и их изменение пользователями осуществляется крайне
редко, что не всегда благоприятным образом сказывается на времени решения
задач. Особого внимания заслуживает открытый свободно распространяемый
пакет OpenFOAM [2], представляющий собой платформу для решения задач
механики сплошной среды. В нем установки «по умолчанию» отсутствуют, и
требуется явно указывать методы решения линейных систем, возникающих в
результате аппроксимации соответствующих уравнений. От того, насколько
удачно будет выбран итерационный метод и, если необходимо, его параметры,
во многом будет зависеть время выполнения расчета.
В пакете OpenFOAM реализованы следующие методы решения СЛАУ:
1) релаксационные методы:
а) метод Гаусса-Зейделя (GS);
б) метод неполного разложения Холецкого (DIC) /
метод неполного LU-разложения (DILU);
в) метод Гаусса-Зейделя с DIC-предобуславливанием (GS + DIC) /
метод Гаусса-Зейделя с DILU-предобуславливанием (GS + DILU);
2) проекционные методы:
а) метод сопряжённых градиентов (CG) / метод бисопряжённых градиентов
(BiCG);
б) метод бисопряжённых градиентов со стабилизацией (BiCGStab) [3];
в) метод обобщенных минимальных невязок (GMRES);
3) многосеточные методы: геометро-алгебраический многосеточный метод
1
Работа проводилась при финансовой поддержке Министерства образования
и науки Российской Федерации
71
(GAMG) и его модификации [4];
72
4) диагональный решатель (для диагональных систем).
Запись вида A / B означает, что метод A является частным случаем метода B
для систем с симметричной матрицей (соответственно для систем с
несимметричной матрицей применяется метод B). Курсивом выделены
группы методов, которые могут быть использованы как предобуславливатели.
Отметим, что методы BiCGStab и GMRES реализованы только в расширенной
версии пакета (OpenFOAM-ext [5]).
Часто выбору итерационного метода решения СЛАУ уделяют мало внимания,
оставляя неизменными те настройки, которые в соответствующих
программных комплексах установлены по умолчанию, либо используются в
похожих примерах в руководствах, как это часто делается при работе с
OpenFOAM. Однако даже при решении простых тестовых задач время счёта
можно сократить в несколько раз, если выбирать методы с учётом специфики
конкретной задачи.
Целью данной работы является построение методики анализа вычислительной
эффективности итерационных методов решения СЛАУ, возникающих при
численном решении разностных аналогов уравнений механики сплошной
среды, а также методики априорного выбора метода решения таких СЛАУ,
обладающего высокой вычислительной эффективностью.
С практической точки зрения наиболее важной характеристикой
итерационного метода решения СЛАУ является время, затрачиваемое на
решение системы с заданной точностью. В то же время при теоретическом
исследовании итерационных методов основное внимание, как правило,
уделяется скорости сходимости, т.е. скорости уменьшения ошибки. Однако
опыт применения различных итерационных методов при решении задач
механики сплошной среды показывает, что прямой связи между скоростью
сходимости метода и его вычислительной эффективностью для конкретной
задачи может не быть. При этом речь идет не только о различающихся
трудоемкостях выполнения одной итерации при использовании различных
методов, но и о необходимости учета таких характеристик, как коэффициент
сглаживания и показатель «затратности», которые определяются не только
численным методом, но и решаемой задачей.
2. Постановка задачи
Рассмотрим систему линейных алгебраических уравнений
Ax  b ( x, b  R N d , det A  0 ),
возникающую
при
дифференциальный
дискретизации
либо
уравнения
(1)
Lu  f , где
интегро-дифференциальный
оператор,
L
f
–
является разреженной и не обладает специальными свойствами (симметрией,
положительной определенностью).
В данной работе в качестве примера будет рассмотрена модельная задача о
расчёте течения вязкой несжимаемой жидкости в квадратной каверне, которое
описывается уравнениями неразрывности и Навье-Стокса:

  v  0,
 


 v
1 
 t  (v  )v  p  Re v .

Здесь Re – число Рейнольдса, v – скорость течения, p – давление; все
величины являются безразмерными.
Задача решается в пространственной области ( x; y )  [0;1]  [0;1] при t  0 .
На границах каверны (при x  0 , x  1 , y  0 и y  1 ) выполняются
граничные условия прилипания и равенства нулю нормальной производной
давления, при этом верхняя «крышка» каверны ( y  1 ) движется с

постоянной скоростью v  (1;0) , остальные «стенки» неподвижны.
Давление в такой постановке определяется с точностью до константы,
поэтому для выделения единственного решения в нижнем правом углу
каверны давление полагается постоянным и равным константе. Разностная
схема для решения задачи строится интегро-интерполяционным методом LSSTAG [6] на прямоугольной структурированной сетке.
На каждом шаге расчёта по времени решение задачи сводится к решению
уравнения Гельмгольца
(   k 2 ) u1  f1
(2)
для прогноза скорости u1 и уравнения Пуассона
u 2  f 2
(3)
для поправки давления u 2 (здесь  – оператор Лапласа). Разностные аналоги
уравнений (2) и (3) представляют собой системы линейных алгебраических
уравнений вида (1). Как показывает практика, при решении этих СЛАУ
итерационными методами вычислительная сложность процедуры решения
разностного аналога уравнения Пуассона (3) оказывается существенно выше,
чем для уравнения Гельмгольца (2).
–
известная функция, u – искомая функция. Будем считать, что матрица A
73
T
74
3. Предобуславливание
Поскольку любой итерационный метод решения системы линейных
алгебраических уравнений (1) может быть представлен как метод, основанный
на расщеплении, запишем его в виде
Mx n 1  Nx n  b , M  N  A .
n
Здесь x – n -е итерационное приближение к искомому решению x . Пусть
r n  b  Ax n – вектор невязки, z n  x  x n – итерационная ошибка,
p n  x n1  x n – вектор коррекции. Легко заметить, что закон изменения
итерационной ошибки при переходе от итерации к итерации имеет
n 1
 M -1 Nz n . Это соотношение показывает, что скорость
-1
сходимости определяется нормой матрицы M N перехода от итерации к
следующий вид: z
итерации, которую можно оценить сверху при помощи спектрального радиуса
[1, 4]. Ясно, что чем точнее M приближает A , тем выше скорость
сходимости. Таким образом вводится понятие предобуславливания: если M
близка к A в смысле некоторой нормы, она является предобуславливателем.
Также отметим, что любой итерационный метод, основанный на расщеплении,
может быть переписан как метод простой итерации (с параметром   1 ) для
решения предобусловленной системы. Понятно, что можно использовать не
один предобуславливатель, т.е. можно заменить метод простой итерации на
другой итерационный метод и использовать его для решения
предобусловленной системы. Например, на рис. 1 показано убывание нормы
невязки при решении разностного аналога уравнения Пуассона (3) на сетке
16  16 методом Якоби без предобуславливания, методом простой итерации с
ILU-предобуславливанием и методом Якоби с ILU-предобуславливанием [1].
Точность 10
16
Рис. 1. Зависимость нормы невязки от номера итерации
В дальнейших примерах в качестве базового метода будем использовать метод
BiCGStab [3], который относится к методам крыловского типа и обладает
наиболее быстрой и гладкой сходимостью из всех методов данного класса.
Алгоритм метода BiCGStab имеет вид:
r 0 = b  Ax 0 , p 0 = r 0 , r*0 : (r*0 , r 0 )  0
for n = 0,1,, while (|| r n ||2   ), do
n =
выбрана исключительно для наглядности.
(r*0 , r n )
( Ap n , r*0 )
s n = r n   n Ap n
n =
( As n , s n )
( As n , As n )
x n1 = x n   n p n   n s n
r n1 = s n   n As n
n =
 n (r n1 , r*0 )
 n (r n , r*0 )
p n1 = r n1   n ( p n   n Ap n ).
75
76
n
k 1
n
здесь p – вектор коррекции, r – вектор невязки на n -й итерации. При
использовании правого предобуславливания [1] вместо исходной системы
Ax = b последовательно решаются системы AM 1 y = b и Mx = y . Такой
последовательности действий соответствует алгоритм метода BiCGStab с
предобуславливанием:
формуле
u kj = a kj  (1   )l kj u ij [7]. Использование такой модификации
i =1
целесообразно, если имеется серия однотипных задач: на одной из них
подбирается оптимальный параметр  , что позволяет сократить
вычислительную трудоемкость процедуры решения по сравнению с
«базовым» вариантом (при   1 ) примерно на 10 %.
r 0 = b  Ax 0 , p 0 = r 0 , r*0 : (r*0 , r 0 )  0
for n = 0,1,, while (|| r n1 ||2   ), do
y n = M 1 p n
n =
(r*0 , r n )
( Ay n , r*0 )
s n = r n   n Ay n
z n = M 1 s n
n =
( Az n , s n )
( Az n , Az n )
x n1 = x n   n y n   n z n
r n1 = s n   n Az n
n =
 n (r n1 , r*0 )
 n (r n , r*0 )
p n1 = r n1   n ( p n   n Ay n ).
Видно, что отличие от исходного алгоритма состоит только в двух шагах (
y n = M 1 p n и z n = M 1 s n ); при этом явного обращения матрицы M при
проведении вычислений не производится, вместо этого решаются системы
My n  p n и Mz n = s n . На рис. 2 показано, что предобуславливание
позволяет многократно увеличить скорость сходимости метода BiCGStab как
для СЛАУ с симметричной матрицей, так и с несимметричной матрицей.
При этом существенный выигрыш получается не только по числу итераций
(табл. 1), но и по числу умножений (табл. 2). Отметим, что  ILUпредобуславливанием
здесь
называется
модификация
ILUпредобуславливания, для которой элементы матрицы U вычисляются по
77
Рис. 2. Убывание норм невязок при решении различными методами
тестовых СЛАУ
78
809  809 : а) симметричной; б) несимметричной.
Табл. 1. Количество итераций при при решении нескольких тестовых СЛАУ
размерности 809  809 различными итерационными методами.
Метод
Метод Зейделя
Метод верхней
«Тест 1»
3313
Количество итераций
«Тест 2»
«Тест 3»
3069
2987
n
2241
2168
2219
BiCGStab
139
41
138
38
240
30
283
40
40
42
30
42
37 (  =0.15 )
36 (  =0.11 )
27 (  =0.3 )
37 (  =0.1 )
37 (  =0.15 )
36 (  =0.11 )
27 (  =0.3 )
38 (  =0.05 )
BiCGStab+ ILU
BiCGStab+ ILU
(ван дер Ворст)
 u ( x)  u ( x) , x  (0,1) ; u (0)  u (1)  0
Табл. 2. Количество умножений при решении тестовых СЛАУ размерности
809  809 различными итерационными методами ( M e – число ненулевых
элементов в матрице СЛАУ).
Метод
в 1 итерацию
Количество умножений
«Тест 1»
«Тест 2»
«Тест 3»
«Тест 4»
Me
17 754 367
16 152 147
14 388 379
15 038 820
Me
12 775 856
11 794 383
10 443 256
11 068 372
BiCGStab
2M e  11N d
2 726 763
2 680 650
4 447 920
5 341 625
BiCGStab+ ILU
4M e  9 N d
1 177 397
1 076 654
796 470
1 089 320
BiCGStab+ ILU
(ван дер Ворст)
BiCGStab+
5M e  7 N d
1 298 320
1 343 076
892 440
1 285 326
4M e  9 N d
1 062 529
1 019 988
716 823
1 007 621
5M e  7 N d
1 200 946
1 151 208
803 196
1 162 914
Метод Зейделя
Метод верхней
релаксации (
w=1.4 )
ILU
BiCGStab+
ILU
(ван дер Ворст)
и
ошибки z в виде разложения по гармоникам, возникающим из решения
разностного аналога одномерной спектральной задачи
2384
BiCGStab+ ILU
(ван дер Ворст)
ошибки
Для дальнейшего анализа эффективности методов необходимо ввести
некоторые дополнительные понятия. Представим вектор итерационной
«Тест 4»
3015
релаксации ( w=1.4 )
BiCGStab+ ILU
4. Анализ изменения итерационной
коэффициенты усиления гармоник
(4)
на равномерной сетке с шагом h . Собственные векторы матрицы СЛАУ,
соответствующей разностному аналогу задачи (4), – гармоники – имеют вид
 kj  sin kjh ,
j  1, N d , k  1, N d .
k  N d неразрешимы на такой сетке,
поскольку длина волны для любой такой гармоники меньше 2h . Остальные
гармоники (с номерами 1  k  N d ) естественным образом делятся на
низкочастотные
( 1  k  [ N d / 2] )
и
высокочастотные
(
[ N d / 2]  1  k  N d ; запись [] означает целую часть числа).
Таким образом, гармоники с номерами
Низкочастотные гармоники являются длинноволновыми (длина волны больше
4h ), а высокочастотные – коротковолновыми.
Наличие длинноволновых и коротковолновых компонент определяет
жёсткость задачи, заключающуюся в различной скорости убывания
n 1
соответствующих компонент ошибки ( z
 M Nz ) при выполнении
одной итерации. Скорость убывания p-й компоненты ошибки можно оценить
коэффициентом усиления соответствующей гармоники
Nd
g p   Mk
k 1
где
Mk
Mk
1
N
1
N
1
N
k ,

p 1
– k -е собственное число матрицы
-1
n
p  1, N d ,
M 1 N (предполагаем, что все
 k – дискретное преобразование Фурье (ДПФ) [8]
действительные), 
k -го собственного вектора матрицы M 1 N , дополненного нулями.
Такие итерационные методы, как метод Якоби, метод Гаусса – Зейделя и т.п.,
редко используются для решения жёстких задач, поскольку за первые
79
80
несколько итераций они практически полностью устраняют высокочастотные
компоненты ошибки, тогда как низкочастотные подавляются очень медленно.
Таким образом, указанные методы действуют как фильтры низких частот,
поэтому их называют сглаживателями, или релаксационными методами, а
максимальный коэффициент усиления высокочастотных гармоник –
коэффициентом сглаживания  sm [4].
Табл. 3. Характеристики метода BiCGStab с ILU-предобуславливанием при
решении разностных аналогов уравнений Гельмгольца и Пуассона.
Разностный
аналог
Размерность
ДПФ превращается в громоздкую задачу, поэтому для рассматриваемой
СЛАУ и итерационного метода важно иметь характеристику, которая
определяется видом дискретизируемого оператора и мало зависит от
величины шага по пространству. Как показывают вычислительные
эксперименты, такой характеристикой является отношение наибольшего
коэффициента усиления низкочастотных гармоник к коэффициенту
сглаживания, который обозначим  и назовём показателем «затратности»
итерационного метода для рассматриваемой задачи, поскольку он показывает,
во сколько раз медленнее подавляются низкочастотные гармоники по
сравнению с высокочастотными, т.е. насколько плохо метод преодолевает
жёсткость задачи. Чем ближе  к единице и чем ближе к нулю спектральный
СЛАУ ( N d )
1
радиус матрицы M N перехода от итерации к итерации, тем лучше метод
подходит для решения данной задачи.
5. Вычислительные эксперименты
Как было отмечено выше, в качестве примеров СЛАУ рассмотрим разностные
аналоги уравнений Гельмгольца (2) и Пуассона (3), полученные при решении
хорошо известной тестовой задачи о моделировании течения в каверне
методом погруженных границ с функциями уровня (LS-STAG) [6]. Как видно
из табл. 3, для решения разностного аналога уравнения Гельмгольца
рассматриваемый метод BiCGStab с ILU-предобуславливанием весьма
эффективен, в то время как для плохо обусловленной системы, полученной
 ( M 1 N )
«затратности» 
при дискретизации уравнения Пуассона, спектральный радиус
матрицы перехода от итерации к итерации и показатель
далеки от нуля и единицы соотвественно, поэтому для численного решения
второго уравнения необходимо подобрать более эффективный метод.
Пуассона
уравнения
N d вычисление собственных чисел, собственных векторов и их
При больших
Гельмгольца
12
56
240
16
64
256
1.00
1.00
1.00
85.88
589.89
3292.74
 ( M 1 N )
1.5  10 9
2.8  10 8
1.7  10 6
0.82
0.97
0.99
 sm
1.7  10 9
7.3  10 8
1.5  10 5
0.13
0.24
0.30

1.01
1.14
1.09
3.34
3.25
3.32
Число
обусловленности
Поскольку наиболее трудоемкой оказывается «борьба» с низкочастотными
составляющими
ошибки,
весьма
перспективной
представляется
многосеточная стратегия. Сначала выполняется несколько сглаживающих
итераций, подавляющих высокочастотные составляющие, а затем
осуществляется переход на более грубую сетку, на которой высокочастотные
гармоники неразрешимы, а половина низкочастотных становятся
высокочастотными для нового сеточного уровня. На самой грубой сетке
задача, как правило, решается методом Гаусса, а затем решение переносится
на исходную сетку. При переносе на подробную сетку выполняются
постсглаживающие итерации.
Многосеточный метод [4, 9–11] является хорошим предобуславливателем.
Здесь для примера рассмотрим модификацию многосеточного метода,
описанную в [7]. В ней присутствуют три параметра, значения которых можно
варьировать. Это параметр релаксации  в используемом в качестве
сглаживателя методе ADLJ, число предсглаживаний n pre и число
постсглаживаний
n post . В [12] указано, что при решении СЛАУ, возникающих
при дискретизации эллиптических уравнений в двумерном случае, для
подобного решателя оптимальными будут следующие параметры:   0.71 ,
n pre  0 , n post  1 или n post  2 . Однако, рассчитав такие характеристики
данного итерационного метода при решении разностного аналога уравнения
Пуассона, как спектральный радиус
81
82
 ( M 1 N )
и показатель «затратности»
,
был сделан вывод о том, что в двумерном случае наилучшим выбором
являются такие значения параметров:
  0.86 , n pre  1 , n post  1 .
При
них наблюдается наименьший спектральный радиус матрицы перехода от
1
итерации к итерации M N при близком к единице показателе
«затратности». Собственные числа матрицы перехода от итерации к итерации
и коэффициенты усиления гармоник для данных параметров представлены на
рис. 3.
Рис. 4. Временные затраты при решении тестовой задачи о моделировании течения в
каверне методом BiCGStab с различными предобуславливателями
(OpenFOAM,
N d  10000 ).
6. Заключение
Определение коэффициентов усиления гармоник позволяет получить ряд
важных характеристик итерационного метода, по которым можно судить о его
эффективности при решении той или иной задачи. Для определения этих
коэффициентов предлагается использовать дискретное преобразование Фурье.
Как показывают вычислительные эксперименты, для итерационного метода
отношение наибольшего коэффициента усиления низкочастотной гармоники к
наибольшему коэффициенту усиления высокочастотной гармоники, названное
показателем «затратности», определяется видом дискретизируемого оператора
и мало зависит от N d . Благодаря введению такой характеристики, которая
Рис. 3. Собственные числа матрицы перехода от итерации к итерации и
коэффициенты усиления гармоник при решении разностного аналога уравнения
Пуассона на сетке
16  16
предобуславливанием ( 
Отметим, что при
методом BiCGStab с многосеточным
 0.86 , n pre  1 , n post  1 ).
N d  17760 и ограничении на норму невязки   10 6 на
первом шаге по времени система решается методом BiCGStab с ILUпредобуславливанием за 158 итераций, а методом BiCGStab с многосеточным
прелобуславливанием с указанными выше параметрами – за 9. Для сравнения,
если использовать параметры, приведенные в [12] (   0.71 , n pre  0 ,
n post  1 ), число итераций увеличивается до 11, при иных параметрах (
  0.10 , n pre  0 , n post  2 ), число итераций увеличивается до 23.
может быть вычислена при решении задачи малой размерности (на грубой
сетке), стало возможным определение оптимальных параметров при
использовании многосеточного предобуславливания.
Работа проводится при финансовой поддержке Министерства образования и
науки Российской Федерации
Список литературы
[1] Saad Y. Iterative Methods for Sparse Linear Systems. N.Y.: PWS Publ., 1996. 547 p.
[2] OpenFOAM // The OpenFOAM Foundation. URL: http://www.openfoam.org (дата
обращения: 30.04.2013).
[3] Van der Vorst H.A. Bi-CGSTAB: a fast and smoothly converging variant of Bi-CG
for solution of non-symmetric linear systems// SIAM J. Sci. Stat. Comp. 1992. № 2.
P. 631-644.
[4] Wesseling P. An introduction to multigrid methods. Chichester: John Willey & Sons
Ltd., 1991. 284 p.
[5] The OpenFOAM® Extend Project. URL: http://www.extend-project.de (дата
обращения: 30.04.2013).
Решая тестовую задачу о моделировании течения в квадратной каверне в
пакете OpenFOAM на сетке 100  100 , наименьшее время счёта также
получается при использовании ILU-предобуславливания для разностного
аналога уравнения Гельмгольца и многосеточного предобуславливания для
разностного аналога уравнения Пуассона (рис. 4).
83
84
[6] Cheny Y., Botella O. The LS-STAG method: A new immersed boundary/level-set
method for the computation of incompressible viscous flows in complex moving
geometries with good conservation properties // J. Comput. Phys. 2010. № 229. P.
1043-1076.
[7] Пузикова В.В. Решение систем линейных алгебраических уравнений методом
BiCGStab с предобуславливанием // Вестник МГТУ им. Н.Э. Баумана.
Естественные науки. 2011. Спец. выпуск «Прикладная математика». С. 124133.
[8] Сергиенко А.Б. Цифровая обработка сигналов. СПб.: Питер, 2002. 608 с.
[9] Федоренко Р.П. Введение в вычислительную физику. М.: Изд-во Моск. физ.техн. ин-та, 1994. 528 с.
[10] Галанин М.П., Савенков Е.Б. Методы численного анализа математических
моделей. М.: Изд-во МГТУ им. Н.Э. Баумана, 2010. 590 с.
[11] Ольшанский М.А. Лекции и упражнения по многосеточным методам. М.: Издво ЦПИ при механико-математическом факультете МГУ им. М.В. Ломоносова,
2003. 163 с.
[12] Van Kan J., Vuik C., Wesseling P. Fast pressure calculation for 2D and 3D time
dependent incompressible flow // Numer. Linear Algebra Appl. 2000. № 7. P. 429447.
85
1/--страниц
Пожаловаться на содержимое документа