Решение систем линейных уравнений в симуляторах

Қысқаша негізгісі

Каждый неявный шаг — система $\mathbf{A}\mathbf{x} = \mathbf{b}$ , по уравнению на ячейку (или на каждую фазу ячейки). В строке — диагональ и по числу на каждую связь: семь чисел в трёхмерной декартовой сетке, остальное — нули.
Матрица уравнения давления симметрична и имеет диагональное преобладание благодаря ёмкостному члену $S/\Delta t$ . Равномерный сдвиг давления — её собственный вектор с наименьшим собственным числом $S/\Delta t$ , поэтому с ростом шага обусловленность ухудшается. Якобиан полностью неявной схемы несимметричен и состоит из блоков.
Метод Гаусса исключает неизвестные и решает треугольную систему обратной подстановкой. Исключение заполняет ленту матрицы: работа порядка $nb^2$ , память — $nb$ ; в трёхмерной сетке ширина ленты — целый слой ячеек, и для больших моделей ленточные методы непригодны.
Якоби, Зейдель и верхняя релаксация уточняют каждое неизвестное по его уравнению. Они сходятся, если спектральный радиус матрицы перехода меньше единицы, но на больших сетках медленно: Якоби и Зейделю нужно $\propto N^2$ итераций, верхней релаксации — $\propto N$ .
Методы Крылова (сопряжённые градиенты, GMRES) строят лучшее приближение из векторов $\mathbf{b}$ , $\mathbf{A}\mathbf{b}$ , $\mathbf{A}^2\mathbf{b}$ , …; у сопряжённых градиентов число итераций растёт как $\sqrt\kappa$ . Предобусловливатель — грубое, но дешёвое подобие $\mathbf{A}$ — уменьшает его в разы, особенно в неоднородных моделях.
Итерации останавливают по относительной невязке. Невязка уравнения давления — небаланс ячейки в м³/сут: недорешённая система «создаёт» флюид, и её погрешность прямо входит в материальный баланс.

Негізгі қатынастар

\left(\frac{S_i}{\Delta t} + \sum_j a_{ij}\right)\delta_i - \sum_j a_{ij}\,\delta_j = -q_i

строка системы давления неявного шага (давления в начале шага одинаковы)

W_{\text{лента}} \approx n\,b^2,\qquad W_{\text{итер}} \approx \nu\,c\,n

работа ленточного исключения и итерационного метода

\nu \approx \frac{\ln(1/\varepsilon)}{-\ln\rho(\mathbf{G})}

итераций, чтобы уменьшить ошибку в $1/\varepsilon$ раз

\frac{\lVert\mathbf{x} - \mathbf{x}^{(\nu)}\rVert}{\lVert\mathbf{x}\rVert} \le \kappa\,\frac{\lVert\mathbf{r}^{(\nu)}\rVert}{\lVert\mathbf{b}\rVert}

малая невязка — ещё не малая ошибка, если $\kappa$ велико

Жиі қателіктер

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

Сабақтан кейін сіз

описать структуру матрицы гидродинамической модели: разреженность, ленту, симметрию, диагональное преобладание и зависимость обусловленности от шага по времени;
решить небольшую систему методом Гаусса, объяснить LU-разложение, заполнение и роль нумерации ячеек и оценить работу ленточного метода;
выполнить итерации методов Якоби, Гаусса — Зейделя и последовательной верхней релаксации и объяснить, от чего зависит их сходимость;
объяснить идею методов подпространств Крылова и предобусловливания — неполного LU-разложения, двухступенчатого CPR, многосеточных методов — и почему они нужны неоднородным моделям;
выбрать критерий остановки, связать невязку с материальным балансом и объяснить, почему решатель определяет время расчёта.

Откуда берутся системы уравнений и как они устроены

Неявная схема (явная и неявная схемы) превращает каждый шаг по времени в систему уравнений: перетоки считают по давлениям конца шага, ещё неизвестным, и уравнения всех ячеек решают совместно. В схеме IMPES это одна линейная система для давлений — по уравнению на ячейку. В полностью неявной схеме нелинейные уравнения решают методом Ньютона, и на каждой итерации — линейную систему с якобианом, по уравнению на каждую фазу ячейки (невязки и якобиан). За расчёт прогноза набираются сотни шагов по несколько итераций — тысячи систем с сотнями тысяч или миллионами неизвестных.

Запишем систему давления одного неявного шага для однофазной слабосжимаемой жидкости. Неизвестные — изменения давления за шаг $\delta_i = p_i^{n+1} - p_i^n$ :

Строка системы давления: уравнение ячейки на неявном шаге

\left(\frac{S_i}{\Delta t} + \sum_j a_{ij}\right)\delta_i - \sum_j a_{ij}\,\delta_j = -q_i + \sum_j a_{ij}\left(p_j^{n} - p_i^{n}\right)

мұндағы:

$\delta_i,\ \delta_j$ изменения давления ячейки и её соседей за шаг — неизвестные, МПа
$S_i$ упругая ёмкость ячейки $V_{p,i}c_t$ , м³/МПа
$a_{ij}$ коэффициент связи $0{,}0853\,T_{ij}/\mu$ , м³/(сут·МПа)
$q_i$ отбор скважиной в пластовых условиях (закачка — со знаком минус), м³/сут
$p_i^{n},\ p_j^{n}$ давления в начале шага, МПа
$\Delta t$ шаг по времени, сут

Бірліктер: метрические промысловые

Мына жағдайда орындалады

линейное уравнение давления: свойства флюида и породы на шаге постоянны (нелинейный случай решают итерациями, и на каждой — такая же система);
сумма — по всем связям ячейки: соседям по граням и несоседним соединениям; скважина с заданным забойным давлением добавляет к диагонали свой коэффициент.

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

Разреженность (разреженная матрица). В строке столько ненулевых чисел, сколько у ячейки связей, плюс диагональ: три в одномерной полосе, пять в двумерной и семь в трёхмерной декартовой сетке; несоседние соединения и скважины добавляют немногие.
Лента. При естественной нумерации — сначала вдоль оси $x$ , затем по строкам и по слоям — соседи ячейки отстоят от неё по номеру на 1, на $N_x$ и на $N_xN_y$ . Ненулевые числа лежат на семи диагоналях, а ширина ленты — целый слой ячеек.
Симметрия. Коэффициент связи одинаков для пар $(i, j)$ и $(j, i)$ , поэтому матрица давления симметрична. Якобиан полностью неявной схемы несимметричен — подвижность берут из ячейки выше по потоку (из чего складывается якобиан), — и каждой ячейке в нём отвечает блок 2 × 2 или 3 × 3.
Диагональное преобладание. Диагональ больше суммы модулей остальных чисел строки на $S_i/\Delta t$ (и на коэффициент скважины с заданным забойным давлением). Симметричная матрица с таким преобладанием и положительной диагональю положительно определена: решение единственно, а многие методы гарантированно сходятся.

Сумма чисел каждой строки равна ёмкостному члену $S_i/\Delta t$ : коэффициент связи входит в неё дважды — с плюсом на диагонали и с минусом вне её. Поэтому при одинаковой ёмкости ячеек равномерное изменение давления во всём пласте — собственный вектор матрицы:

Равномерный сдвиг давления — собственный вектор матрицы

\mathbf{A} = \frac{S}{\Delta t}\,\mathbf{I} + \mathbf{C},\qquad \mathbf{C}\,\mathbf{1} = 0,\qquad \lambda_{\min}(\mathbf{A}) = \frac{S}{\Delta t},\qquad \kappa = 1 + \frac{\Delta t}{S}\,\lambda_{\max}(\mathbf{C})

мұндағы:

$\mathbf{A}$ матрица системы давления, м³/(сут·МПа)
$\mathbf{C}$ матрица связей: на диагонали $\sum_j a_{ij}$ , вне её $-a_{ij}$ , м³/(сут·МПа)
$\mathbf{I},\ \mathbf{1}$ единичная матрица и вектор из единиц — одинаковое изменение давления во всех ячейках, —
$S$ упругая ёмкость ячейки (одинаковая у всех ячеек), м³/МПа
$\lambda_{\min},\ \lambda_{\max}$ наименьшее и наибольшее собственные числа, м³/(сут·МПа)
$\kappa$ число обусловленности $\lambda_{\max}(\mathbf{A})/\lambda_{\min}(\mathbf{A})$ , —

Бірліктер: метрические промысловые; $\kappa$ безразмерно

Мына жағдайда орындалады

одинаковая ёмкость ячеек; внешние границы закрыты, скважины работают с заданным дебитом;
скважина с заданным забойным давлением или граница с постоянным давлением добавляет к диагонали и увеличивает $\lambda_{\min}$ .

Связи не «видят» равномерного сдвига давления — его определяет только ёмкость. Чем крупнее шаг, тем меньше $S/\Delta t$ и тем хуже обусловлена система; в пределе несжимаемой жидкости с заданными дебитами она вырождена — давление определено с точностью до постоянной. Именно такие гладкие, общие для большой части модели составляющие ошибки труднее всего убрать простыми итерациями.

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

\delta_1

на местах (2, 3) и (3, 2) появляются новые ненулевые элементы — заполнение.

Есептеу мысалыСистема уравнений для участка из четырёх ячеек

Бастапқы деректер (шартты)

Участок: 2 × 2 ячейки 300 × 300 × 12 м: 1 — юго-западная, 2 — юго-восточная, 3 — северо-западная, 4 — северо-восточная; внешние грани закрыты
Проницаемость ячеек 1–4: 120; 40; 200; 80 мД
Порода и жидкость: $\phi$ = 0,20; $c_t = 1{,}5\cdot10^{-3}$ МПа⁻¹; вязкость 1,0 мПа·с
Скважины и шаг: закачка 400 м³/сут в ячейку 1, отбор 400 м³/сут из ячейки 4 (пластовые условия); неявный шаг $\Delta t$ = 10 сут от одинакового давления

Шешім

Проводимости — две половины ячеек последовательно (проводимость связи); у одинаковых ячеек с гранью $A = 300\cdot12 = 3600$ м² и половиной ячейки 150 м $T = 3600/(150/k_i + 150/k_j)$ : $T_{12} = 720$ , $T_{34} = 1371{,}43$ , $T_{13} = 1800$ , $T_{24} = 640$ мД·м.
Коэффициенты связи $a = 0{,}0853\,T/1{,}0$ : $a_{12} = 61{,}416$ , $a_{34} = 116{,}983$ , $a_{13} = 153{,}540$ , $a_{24} = 54{,}592$ м³/(сут·МПа). Ёмкость ячейки $S = 300\cdot300\cdot12\cdot0{,}20\cdot1{,}5\cdot10^{-3} = 324$ м³/МПа, $S/\Delta t = 32{,}4$ .
Диагональ — $S/\Delta t$ плюс коэффициенты связей ячейки: 247,356; 148,408; 302,923; 203,975.
Уравнения ячеек (давления в начале шага одинаковы, поэтому справа только скважины): 1: $247{,}356\,\delta_1 - 61{,}416\,\delta_2 - 153{,}540\,\delta_3 = 400$ ; 2: $-61{,}416\,\delta_1 + 148{,}408\,\delta_2 - 54{,}592\,\delta_4 = 0$ ; 3: $-153{,}540\,\delta_1 + 302{,}923\,\delta_3 - 116{,}983\,\delta_4 = 0$ ; 4: $-54{,}592\,\delta_2 - 116{,}983\,\delta_3 + 203{,}975\,\delta_4 = -400$ .
Собственные числа матрицы: 32,4; 176,4; 221,0; 472,9 — наименьшее равно $S/\Delta t$ ; число обусловленности $472{,}9/32{,}4 = 14{,}6$ .

Жауап: Матрица 4 × 4 симметрична; в каждой строке диагональ больше суммы модулей остальных чисел на 32,4; четыре элемента из шестнадцати — нули: ячейки 1 и 4, 2 и 3 гранью не связаны.

Бұл нені білдіреді. Нули на местах (1, 4) и (2, 3) — частный случай главного свойства: в строке только свои связи. В модели из миллиона ячеек в строке те же семь чисел, а нулей — почти миллион. Наименьшее собственное число равно $S/\Delta t$ — это равномерный сдвиг давления. При шаге 100 сут оно стало бы вдесятеро меньше, а число обусловленности выросло бы почти вдесятеро (таблица ниже). Данные примера условные.

Прямые методы: исключение Гаусса и LU-разложение

Метод Гаусса (Gaussian elimination) решает систему в два хода. Прямой ход исключает неизвестные по очереди: из всех уравнений ниже первого вычитают первое, умноженное на такие числа, чтобы исчезло первое неизвестное; затем так же второе — и так до последнего. Матрица становится верхнетреугольной. Обратный ход — подстановка снизу вверх: последнее уравнение содержит одно неизвестное, каждое предыдущее — одно новое.

Шаг прямого хода метода Гаусса

m_{ki} = \frac{A_{ki}}{A_{ii}},\qquad A_{kl} \leftarrow A_{kl} - m_{ki}\,A_{il},\qquad b_k \leftarrow b_k - m_{ki}\,b_i\qquad (k > i)

мұндағы:

$A_{kl},\ b_k$ элементы матрицы и правой части, которые меняются при исключении, в единицах системы
$i$ номер исключаемого неизвестного (ведущая строка), —
$m_{ki}$ множитель строки $k$ ; множители образуют нижнетреугольную матрицу $\mathbf{L}$ разложения, —

Бірліктер: в единицах системы

Мына жағдайда орындалады

ведущий элемент $A_{ii}$ не равен нулю; при диагональном преобладании это гарантировано, и перестановки строк не нужны;
строки, в которых $A_{ki} = 0$ , не меняются.

После исключения всех неизвестных матрица становится верхнетреугольной ( $\mathbf{U}$ ), а множители образуют нижнетреугольную матрицу $\mathbf{L}$ с единицами на диагонали: $\mathbf{A} = \mathbf{L}\mathbf{U}$ — LU-разложение. Его делают один раз, а затем для любой правой части решают две треугольные системы. Для симметричной положительно определённой матрицы применяют разложение Холецкого $\mathbf{A} = \mathbf{G}\mathbf{G}^{\mathsf T}$ — вдвое дешевле.

Обратная подстановка

x_n = \frac{\tilde b_n}{U_{nn}},\qquad x_i = \frac{1}{U_{ii}}\left(\tilde b_i - \sum_{l > i} U_{il}\,x_l\right),\quad i = n - 1,\ \ldots,\ 1

мұндағы:

$U_{il}$ элементы верхнетреугольной матрицы после прямого хода, в единицах системы
$\tilde b_i$ правая часть после прямого хода, в единицах системы
$x_i$ неизвестные, в единицах системы

Бірліктер: в единицах системы

Мына жағдайда орындалады

верхнетреугольная система с ненулевой диагональю.

Для трёхдиагональной матрицы оба хода и составляют метод прогонки: число действий пропорционально числу неизвестных.

Есептеу мысалыМетод Гаусса для участка из четырёх ячеек

Бастапқы деректер (шартты)

Система: из предыдущего примера; коэффициенты — с тремя знаками после запятой, множители — с пятью

Шешім

Исключаем $\delta_1$ из строки 2. Множитель $61{,}416/247{,}356 = 0{,}24829$ ; прибавляем к строке 2 строку 1, умноженную на него: диагональ $148{,}408 - 61{,}416\cdot0{,}24829 = 133{,}159$ , правая часть $400\cdot0{,}24829 = 99{,}316$ , а на месте нуля при $\delta_3$ появляется $-153{,}540\cdot0{,}24829 = -38{,}122$ — это заполнение.
Строка 3: множитель $153{,}540/247{,}356 = 0{,}62072$ ; диагональ $302{,}923 - 153{,}540\cdot0{,}62072 = 207{,}618$ , при $\delta_2$ появляется $-61{,}416\cdot0{,}62072 = -38{,}122$ , правая часть $400\cdot0{,}62072 = 248{,}288$ . Строку 4 исключение $\delta_1$ не меняет: в ней $\delta_1$ нет.
Исключаем $\delta_2$ с помощью новой строки 2. Строка 3: множитель $38{,}122/133{,}159 = 0{,}28629$ ; диагональ $207{,}618 - 38{,}122\cdot0{,}28629 = 196{,}704$ , при $\delta_4$ : $-(116{,}983 + 54{,}592\cdot0{,}28629) = -132{,}612$ , правая часть $248{,}288 + 99{,}316\cdot0{,}28629 = 276{,}721$ .
Строка 4: множитель $54{,}592/133{,}159 = 0{,}40998$ ; диагональ $203{,}975 - 54{,}592\cdot0{,}40998 = 181{,}593$ , при $\delta_3$ : $-(116{,}983 + 38{,}122\cdot0{,}40998) = -132{,}612$ , правая часть $-400 + 99{,}316\cdot0{,}40998 = -359{,}282$ .
Исключаем $\delta_3$ из строки 4: множитель $132{,}612/196{,}704 = 0{,}67417$ ; диагональ $181{,}593 - 132{,}612\cdot0{,}67417 = 92{,}190$ , правая часть $-359{,}282 + 276{,}721\cdot0{,}67417 = -172{,}725$ . Матрица стала верхнетреугольной.
Обратная подстановка: $\delta_4 = -172{,}725/92{,}190 = -1{,}87358$ ; $\delta_3 = (276{,}721 + 132{,}612\,\delta_4)/196{,}704 = 0{,}14368$ ; $\delta_2 = (99{,}316 + 38{,}122\,\delta_3 + 54{,}592\,\delta_4)/133{,}159 = 0{,}01886$ ; $\delta_1 = (400 + 61{,}416\,\delta_2 + 153{,}540\,\delta_3)/247{,}356 = 1{,}71097$ МПа.
Проверка: решение без округлений — 1,71099; 0,01887; 0,14370; −1,87356 МПа; округления в цепочке изменили ответ не больше чем на 0,00002 МПа.

Жауап: За шаг давление в ячейке с закачкой выросло на 1,71 МПа, в ячейке с отбором упало на 1,87 МПа, в ячейках 2 и 3 изменилось на 0,019 и 0,14 МПа.

Бұл нені білдіреді. Исключив ячейку 1, мы связали её соседей — ячейки 2 и 3, которые гранью не связаны: в матрице появился новый ненулевой элемент. Так происходит всегда, и в большой сетке заполнение занимает всю ленту (рисунок ниже). Обратите внимание на $\delta_2$ : это небольшая разность больших чисел, и при вычислениях с малым числом знаков она теряла бы точность; машинная арифметика держит около 16 знаков, но плохо обусловленные системы теряют их тем больше, чем больше $\kappa$ . Данные примера условные.

ЕсептеуМатрица уравнения давления сетки 4 × 3 при естественной нумерации (вдоль строк по 4 ячейки) и множители её LU-разложения: ненулевые элементы матрицы лежат на пяти диагональных линиях, а исключение заполняет почти всё между крайними — ленту шириной в строку сетки (36 новых элементов). При нумерации вдоль короткой стороны (по 3 ячейки) лента уже и заполнение меньше — 24 элемента.

Заполнение и нумерация. Исключение неизвестного связывает всех его ещё не исключённых соседей друг с другом, поэтому в сетке при естественной нумерации заполняется почти вся лента между крайними ненулевыми диагоналями. Работа и память растут с её шириной, и порядок нумерации важен: нумеруют сначала вдоль самого короткого направления. Для сетки 10 × 5 нумерация вдоль длинной стороны даёт 648 новых элементов, вдоль короткой — 288. Обратная матрица разреженной системы заполнена почти целиком, поэтому её никогда не строят (якобиан не обращают).

Работа и память ленточного исключения

W \approx n\,b^2,\qquad M \approx n\,b;\qquad b = N_x\ \ (\text{2D}),\qquad b = N_xN_y\ \ (\text{3D})

мұндағы:

$W$ число операций умножения со сложением, —
$M$ число хранимых чисел разложения Холецкого (у LU-разложения — вдвое больше: $\mathbf{L}$ и $\mathbf{U}$ ), —
$n$ число неизвестных, —
$b$ полуширина ленты — наибольшая разность номеров связанных неизвестных, —
$N_x,\ N_y$ число ячеек вдоль осей, по которым идёт нумерация (быстрее всего — вдоль $x$ ), —

Бірліктер: безразмерные счёты

Мына жағдайда орындалады

главный член оценки при $b \ll n$ ; исключение без перестановок строк;
для разложения Холецкого симметричной матрицы работа примерно вдвое меньше; для плотной матрицы работа — $n^3/3$ .

Работа растёт как квадрат ширины ленты, поэтому в трёхмерной модели нумеруют сначала по вертикали — обычно это самое короткое направление. Даже тогда ширина ленты — тысячи, и для моделей в сотни тысяч ячеек ленточное исключение слишком дорого (пример в последнем разделе).

Для разреженных матриц есть и более экономные прямые методы: неизвестные перенумеровывают так, чтобы заполнение было меньше (методы минимальной степени, вложенных сечений — nested dissection), и хранят только ненулевые числа. В двумерных задачах они работают хорошо; в трёхмерных их работа растёт примерно как квадрат числа неизвестных, и для моделей в миллионы ячеек, которые решают тысячи раз за расчёт, они слишком дороги. В симуляторах прямые методы остаются для небольших систем — двумерных разрезов, уравнений скважин, грубых уровней многосеточных методов — и как основа предобусловливателей.

Простые итерационные методы: Якоби, Зейдель, верхняя релаксация

Итерационный метод начинает с приближения — например, с нулевых изменений давления (давлений начала шага) — и уточняет его. Простейшая идея: каждое уравнение решить относительно «своего» неизвестного, подставив значения соседей из прошлого приближения, и повторять, пока изменения не станут малы.

Методы Якоби, Гаусса — Зейделя и последовательной верхней релаксации

\begin{aligned} &\text{Якоби:} && x_i^{(\nu+1)} = \frac{1}{A_{ii}}\Bigl(b_i - \sum_{l \ne i} A_{il}\,x_l^{(\nu)}\Bigr)\\ &\text{Зейдель:} && x_i^{(\nu+1)} = \frac{1}{A_{ii}}\Bigl(b_i - \sum_{l < i} A_{il}\,x_l^{(\nu+1)} - \sum_{l > i} A_{il}\,x_l^{(\nu)}\Bigr)\\ &\text{SOR:} && x_i^{(\nu+1)} = x_i^{(\nu)} + \omega\,\bigl(\hat x_i - x_i^{(\nu)}\bigr)\end{aligned}

мұндағы:

$x_i^{(\nu)}$ приближение неизвестного $i$ на итерации $\nu$ , в единицах неизвестных
$A_{il},\ b_i$ элементы матрицы и правой части, в единицах системы
$\hat x_i$ значение, которое дала бы формула Зейделя, в единицах неизвестных
$\omega$ параметр верхней релаксации, $1 < \omega < 2$ , —

Бірліктер: в единицах системы

Мына жағдайда орындалады

$A_{ii} \ne 0$ ; Зейдель и SOR проходят уравнения по порядку номеров и сразу используют новые значения;
сходимость гарантирована при строгом диагональном преобладании (Якоби и Зейдель), а для симметричной положительно определённой матрицы — у Зейделя и у SOR при $0 < \omega < 2$ .

Якоби пересчитывает все неизвестные по значениям прошлой итерации; Зейдель сразу использует уже уточнённые; последовательная верхняя релаксация (successive over-relaxation, SOR) «перешагивает» за поправку Зейделя в $\omega$ раз. Итерация стоит примерно одного умножения матрицы на вектор.

Есептеу мысалыИтерации для участка из четырёх ячеек

Бастапқы деректер (шартты)

Система: участок 2 × 2 из первого примера, шаг 10 сут
Начальное приближение: $\delta = 0$ — давления начала шага
Критерий остановки: относительная невязка $\lVert\mathbf{b} - \mathbf{A}\boldsymbol\delta\rVert/\lVert\mathbf{b}\rVert < 10^{-6}$

Шешім

Якоби, итерация 1: в уравнениях 1 и 4 соседи пока нулевые, $\delta_1 = 400/247{,}356 = 1{,}6171$ , $\delta_4 = -400/203{,}975 = -1{,}9610$ ; в уравнениях 2 и 3 правые части нулевые, и $\delta_2 = \delta_3 = 0$ .
Якоби, итерация 2: $\delta_2 = (61{,}416\cdot1{,}6171 + 54{,}592\cdot(-1{,}9610))/148{,}408$ , то есть $(99{,}316 -107{,}055)/148{,}408 = -0{,}0521$ ; $\delta_3 = (248{,}290 -229{,}404)/302{,}923 = 0{,}0623$ ; $\delta_1$ и $\delta_4$ не изменились — их соседи на итерации 1 были нулевыми.
Зейдель, итерация 1: $\delta_1 = 1{,}6171$ , как у Якоби; затем новое $\delta_1$ сразу идёт в уравнения 2 и 3: $\delta_2 = 99{,}316/148{,}408 = 0{,}6692$ , $\delta_3 = 248{,}290/302{,}923 = 0{,}8196$ ; в уравнении 4 — новые $\delta_2$ и $\delta_3$ : $\delta_4 = (-400 + 36{,}533 + 95{,}879)/203{,}975 = -1{,}3119$ .
Спектральный радиус матрицы Якоби $\rho_J$ = 0,858; у Зейделя $\rho_J^2$ = 0,736; оптимальный параметр $\omega = 2/(1 + \sqrt{1 - \rho_J^2}) = 1{,}32$ даёт $\rho$ = 0,32.
Итераций до невязки $10^{-6}$ : Якоби — 65, Зейдель — 39, верхняя релаксация с $\omega$ = 1,32 — 15; метод сопряжённых градиентов (следующий раздел) — 3.

Жауап: Якоби нужно 65 итераций, Зейделю — 39, верхней релаксации — 15, сопряжённым градиентам — 3.

Бұл нені білдіреді. На первой итерации Зейдель «перескочил»: $\delta_2$ и $\delta_3$ получили 0,6692 и 0,8196 МПа — намного больше решения, потому что в уравнениях 2 и 3 он уже видел новое $\delta_1$ , но ещё не видел отрицательного $\delta_4$ . По первым итерациям методы не сравнивают: скорость в пределе задаёт спектральный радиус, и Зейдель затем сходится вдвое быстрее Якоби. Сопряжённые градиенты закончили за 3 итерации: для такой маленькой системы метод Крылова точен за число шагов, не большее числа неизвестных. Данные примера условные.

Почему итерации сходятся и как быстро? Ошибка приближения $\mathbf{e}^{(\nu)} = \mathbf{x} - \mathbf{x}^{(\nu)}$ на каждой итерации умножается на одну и ту же матрицу перехода $\mathbf{G}$ (дополнительный раздел). Итерации сходятся при любом начальном приближении, если спектральный радиус $\rho(\mathbf{G})$ — наибольший модуль собственного числа — меньше единицы, и в пределе ошибка убывает в $1/\rho$ раз за итерацию:

Скорость сходимости простых итераций

\nu \approx \frac{\ln(1/\varepsilon)}{-\ln\rho(\mathbf{G})},\qquad \rho_{GS} = \rho_J^2,\qquad \omega_{opt} = \frac{2}{1 + \sqrt{1 - \rho_J^2}},\qquad \rho_{SOR} = \omega_{opt} - 1

мұндағы:

$\nu$ число итераций, чтобы уменьшить ошибку в $1/\varepsilon$ раз, —
$\rho(\mathbf{G})$ спектральный радиус матрицы перехода метода, —
$\rho_J,\ \rho_{GS},\ \rho_{SOR}$ спектральные радиусы методов Якоби, Зейделя и верхней релаксации, —
$\omega_{opt}$ оптимальный параметр верхней релаксации, —

Бірліктер: безразмерные

Мына жағдайда орындалады

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

Зейдель сходится вдвое быстрее Якоби, а верхняя релаксация с оптимальным $\omega$ — в разы быстрее Зейделя. Но $\rho_J$ стремится к единице при измельчении сетки и при росте шага по времени: медленнее всего затухают ошибки, гладкие на масштабе многих ячеек.

Участок 2 × 2: шаг по времени, обусловленность и число итераций

Шаг $\Delta t$ , сут: 1
$\rho_J$: 0,387
$\kappa$: 2,36
Якоби; Зейдель; SOR: 12; 8; 6
Сопряжённые градиенты: 3

Шаг $\Delta t$ , сут: 10
$\rho_J$: 0,858
$\kappa$: 14,6
Якоби; Зейдель; SOR: 65; 39; 15
Сопряжённые градиенты: 3

Шаг $\Delta t$ , сут: 100
$\rho_J$: 0,984
$\kappa$: 137
Якоби; Зейдель; SOR: 464; 287; 46
Сопряжённые градиенты: 3

Шаг $\Delta t$ , сут	$\rho_J$	$\kappa$	Якоби; Зейдель; SOR	Сопряжённые градиенты
1	0,387	2,36	12; 8; 6	3
10	0,858	14,6	65; 39; 15	3
100	0,984	137	464; 287; 46	3

Итерации от нулевого приближения до относительной невязки $10^{-6}$ ; SOR — с оптимальным $\omega$ (1,04; 1,32; 1,69). Величина $\kappa - 1$ растёт строго пропорционально шагу, как предсказывает формула для равномерного сдвига давления, а с ней — и число простых итераций.

На большой сетке простым методам ещё труднее. За одну итерацию Якоби поправка доходит от ячейки только до её соседей, а чтобы исправить гладкую ошибку — например, общий уровень давления в большом блоке пласта, — сведения должны пройти через всю сетку. Поэтому число итераций Якоби и Зейделя растёт как квадрат числа ячеек вдоль стороны, а верхней релаксации — пропорционально ему (таблица в следующем разделе). В симуляторах долго применяли блочные варианты, например линейную верхнюю релаксацию (line SOR, LSOR): целую строку ячеек решают сразу методом прогонки. Сегодня простые итерации работают главным образом внутри более сильных методов — как сглаживатели в многосеточных методах и предобусловливателях.

Методы подпространств Крылова и предобусловливание

Простые итерации используют только последнее приближение. Методы подпространств Крылова (Krylov subspace methods) помнят историю: приближение на итерации $\nu$ ищут как лучшую комбинацию векторов $\mathbf{b}$ , $\mathbf{A}\mathbf{b}$ , …, $\mathbf{A}^{\nu - 1}\mathbf{b}$ , которые получаются умножением матрицы на вектор. Для симметричной положительно определённой матрицы — уравнения давления — это метод сопряжённых градиентов (conjugate gradients, CG): он минимизирует ошибку в «энергетической» норме и в точной арифметике заканчивает работу не больше чем за $n$ итераций. В примере с участком 2 × 2 — за три: в правой части нет составляющей равномерного сдвига (закачка равна отбору), и оставшихся направлений три. Для несимметричных якобианов применяют GMRES и ORTHOMIN — они хранят историю направлений, и её длину ограничивают, отбрасывая старые направления или периодически начиная заново, — или BiCGStab, который обходится короткой историей, но сходится менее гладко.

Сходимость метода сопряжённых градиентов

\lVert\mathbf{e}^{(\nu)}\rVert_{\mathbf{A}} \le 2\left(\frac{\sqrt\kappa - 1}{\sqrt\kappa + 1}\right)^{\nu}\lVert\mathbf{e}^{(0)}\rVert_{\mathbf{A}},\qquad \nu \lesssim \frac{\sqrt\kappa}{2}\,\ln\frac{2}{\varepsilon}

мұндағы:

$\mathbf{e}^{(\nu)}$ ошибка приближения после $\nu$ итераций, в единицах неизвестных
$\lVert\mathbf{e}\rVert_{\mathbf{A}}$ энергетическая норма $\sqrt{\mathbf{e}^{\mathsf T}\mathbf{A}\mathbf{e}}$ , —
$\kappa$ число обусловленности матрицы, —
$\varepsilon$ во сколько раз нужно уменьшить ошибку (обратная величина), —

Бірліктер: безразмерные

Мына жағдайда орындалады

симметричная положительно определённая матрица; точная арифметика;
оценка сверху: на практике сходимость часто быстрее, особенно если собственные числа собраны в группы.

Число итераций растёт как корень из числа обусловленности: для модельной задачи — пропорционально числу ячеек вдоль стороны, как у верхней релаксации, но без подбора параметра. Отсюда смысл предобусловливания — уменьшить $\kappa$ .

Предобусловливание (предобусловливатель, preconditioner). Вместо $\mathbf{A}\mathbf{x} = \mathbf{b}$ решают равносильную систему, в которой матрица ближе к единичной:

Предобусловленная система

\mathbf{M}^{-1}\mathbf{A}\,\mathbf{x} = \mathbf{M}^{-1}\mathbf{b},\qquad \mathbf{M} \approx \mathbf{A},\qquad \kappa\left(\mathbf{M}^{-1}\mathbf{A}\right) \ll \kappa(\mathbf{A})

мұндағы:

$\mathbf{M}$ предобусловливатель: подобие $\mathbf{A}$ , система с которым решается дёшево, в единицах $\mathbf{A}$
$\kappa$ число обусловленности, —

Бірліктер: в единицах системы

Мына жағдайда орындалады

матрицу $\mathbf{M}^{-1}$ не строят: на каждой итерации решают систему $\mathbf{M}\mathbf{z} = \mathbf{r}$ ;
для метода сопряжённых градиентов $\mathbf{M}$ должна быть симметричной положительно определённой.

Хороший предобусловливатель — компромисс: чем ближе $\mathbf{M}$ к $\mathbf{A}$ , тем меньше итераций, но тем дороже каждая. Крайние случаи — $\mathbf{M} = \mathbf{A}$ (одна итерация, но это и есть прямое решение) и диагональ $\mathbf{A}$ (дёшево, но помогает мало).

Неполное LU-разложение (incomplete LU, ILU): исключение Гаусса, в котором заполнение отбрасывают; в ILU(0) множители хранят только на местах ненулевых элементов $\mathbf{A}$ , в ILU(k) — и часть заполнения. Дёшево и особенно полезно для неоднородных моделей.
Вложенная факторизация (nested factorization) — предобусловливатель, придуманный для гидродинамических моделей: он использует вложенный порядок связей по осям сетки (внутри линии ячеек, между линиями, между плоскостями) и широко применяется в симуляторах с декартовыми и угловыми сетками.
Двухступенчатое предобусловливание CPR (constrained pressure residual) для полностью неявной схемы. Сначала из уравнений ячейки составляют уравнение для давления и решают его отдельно: эта часть задачи «эллиптическая» — каждая ячейка чувствует все остальные, и ошибка давления гладкая и общая для всей модели. Затем остаток ошибки, связанный в основном с насыщенностями и местный, убирают дешёвым ILU для всей системы.
Многосеточные методы (multigrid; алгебраические, AMG, строят грубые уровни прямо по матрице): гладкую ошибку, которую простые итерации убирают медленно, исправляют на грубых сетках, где она перестаёт быть гладкой. Число итераций почти не зависит от размера сетки, поэтому AMG часто решает уравнение давления на первой ступени CPR.

Модельная задача: число итераций до относительной невязки

10^{-6}

Метод: Якоби
10 × 10: 258
20 × 20: 986
40 × 40: 3729
20 × 20, неоднородный: 2575

Метод: Зейдель
10 × 10: 130
20 × 20: 494
40 × 40: 1866
20 × 20, неоднородный: 1297

Метод: SOR с оптимальным $\omega$
10 × 10: 28
20 × 20: 53
40 × 40: 102
20 × 20, неоднородный: 92

Метод: Сопряжённые градиенты
10 × 10: 13
20 × 20: 33
40 × 40: 66
20 × 20, неоднородный: 299

Метод: Сопряжённые градиенты + ILU(0)
10 × 10: 10
20 × 20: 16
40 × 40: 30
20 × 20, неоднородный: 24

Метод	10 × 10	20 × 20	40 × 40	20 × 20, неоднородный
Якоби	258	986	3729	2575
Зейдель	130	494	1866	1297
SOR с оптимальным $\omega$	28	53	102	92
Сопряжённые градиенты	13	33	66	299
Сопряжённые градиенты + ILU(0)	10	16	30	24

Квадратный участок пласта, окружённый водоносным горизонтом с постоянным давлением; несжимаемая жидкость, закачка в центре; итерации от нулевого приближения. Неоднородный вариант — логнормальная проницаемость со стандартным отклонением $\ln k$ , равным 2: она меняется от ячейки к ячейке более чем на 4 порядка. Оптимальное $\omega$ для SOR — 1,53; 1,73; 1,85; 1,82. Данные условные.

Таблица показывает три вещи. При удвоении сетки число итераций Якоби и Зейделя растёт почти вчетверо (Якоби: 258 → 986 → 3729), у верхней релаксации и сопряжённых градиентов — примерно вдвое. Неоднородность замедляет простые методы в 1,7–2,6 раза, а сопряжённые градиенты без предобусловливания — примерно в 9 раз: их скорость задаёт обусловленность, которую контрасты проводимостей ухудшают на порядки. Неполное LU-разложение почти снимает эту трудность: 24 итерации против 299. В однородной задаче ILU(0) уменьшает число итераций лишь в 1,3–2,2 раза и не меняет его роста с размером сетки (10 → 16 → 30); рост убирают многосеточные методы.

ЕсептеуОтносительная невязка модельной задачи 20 × 20 (однородный пласт) по итерациям: Якоби, Зейдель, верхняя релаксация (SOR) с оптимальным

\omega

, сопряжённые градиенты (CG) без предобусловливания и с ILU(0). По вертикали — десятичный логарифм невязки: −6 соответствует

10^{-6}

. Якоби и Зейделю до

10^{-6}

нужно около тысячи и пятисот итераций.

Параллельные вычисления

Модели в миллионы ячеек считают на многих процессорах, и решатель должен хорошо делиться между ними. Умножение матрицы на вектор делится легко: каждый процессор отвечает за свою часть сетки и обменивается лишь значениями на её границах. Труднее с предобусловливателем: неполное LU-разложение по природе последовательно, поэтому его строят внутри частей сетки (декомпозиция области), жертвуя частью эффективности; многосеточные методы делятся лучше.

Когда остановить итерации: невязка, ошибка и баланс

Критерий остановки и связь невязки с ошибкой

\frac{\lVert\mathbf{r}^{(\nu)}\rVert}{\lVert\mathbf{b}\rVert} \le \varepsilon,\qquad \mathbf{r}^{(\nu)} = \mathbf{b} - \mathbf{A}\,\mathbf{x}^{(\nu)},\qquad \frac{\lVert\mathbf{x} - \mathbf{x}^{(\nu)}\rVert}{\lVert\mathbf{x}\rVert} \le \kappa\,\frac{\lVert\mathbf{r}^{(\nu)}\rVert}{\lVert\mathbf{b}\rVert}

мұндағы:

$\mathbf{r}^{(\nu)}$ невязка (residual) приближения $\mathbf{x}^{(\nu)}$ , в единицах правой части
$\varepsilon$ допуск — относительная невязка, при которой итерации останавливают, —
$\kappa$ число обусловленности, —
$\lVert\cdot\rVert$ норма вектора (например, корень из суммы квадратов), —

Бірліктер: безразмерные отношения

Мына жағдайда орындалады

невязку можно вычислить, ошибку — нет: точное решение неизвестно; неравенство для ошибки — оценка сверху.

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

У невязки уравнения давления есть физический смысл. Каждое уравнение — баланс ячейки за шаг в кубометрах в сутки, и невязка — объём, «потерянный» или «созданный» в ячейке за сутки. Сумма невязок по всем ячейкам — погрешность материального баланса шага: перетоки между ячейками при сложении взаимно уничтожаются (почему сетка не теряет массу). Поэтому критерии сходимости симуляторов формулируют через нормированные невязки ячеек и общий баланс (критерии сходимости).

Есептеу мысалыЧем грозит грубый допуск: баланс участка из четырёх ячеек

Бастапқы деректер (шартты)

Система: участок 2 × 2, шаг 10 сут; метод Зейделя от нулевого приближения
Допуски: относительная невязка $10^{-2}$ и $10^{-4}$

Шешім

До невязки $10^{-2}$ нужно 9 итераций (невязка 0,0088). Изменения давления: 1,7793; 0,0714; 0,2038; −1,8250 МПа против точных 1,7110; 0,0189; 0,1437; −1,8736 — ошибка до 0,068 МПа.
Ёмкость, умноженная на сумму изменений: $S\sum\delta_i = 324\cdot0{,}2295 = 74$ м³. В точном решении сумма равна нулю — закачка равна отбору. Значит, недорешённые давления «создали» 74 м³ жидкости — 1,9 % закачанных за шаг 4000 м³.
До невязки $10^{-4}$ — 24 итерации: ошибка давления 0,0007 МПа, «созданный» объём 0,8 м³, или 0,02 % закачки.

Жауап: При допуске

10^{-2}

ошибка давлений — до 0,068 МПа, а баланс шага нарушен на 74 м³ (1,9 % закачки); при

10^{-4}

— до 0,0007 МПа и 0,8 м³.

Бұл нені білдіреді. Невязка $10^{-2}$ звучит как «один процент», но это процент от правой части, а не от давлений: изменение давления в ячейке 2, где правая часть нулевая, вышло в 3,8 раза больше точного. Ошибка баланса накапливается от шага к шагу; по ней симулятор и судит, достаточно ли решены уравнения. В полностью неявной схеме допуски линейного решателя согласуют с итерациями Ньютона: на первых итерациях достаточно уменьшить невязку на один-два порядка, на последних — сильнее. Данные примера условные.

Сколько это стоит: почему решатель определяет время расчёта

Работа прямого и итерационного решения

W_{\text{плотн}} \approx \frac{n^3}{3},\qquad W_{\text{лента}} \approx n\,b^2,\qquad W_{\text{итер}} \approx \nu\,c\,n

мұндағы:

$n$ число неизвестных, —
$b$ полуширина ленты, —
$\nu$ число итераций, —
$c$ операций на неизвестное за итерацию: умножение матрицы на вектор, предобусловливатель, скалярные произведения, —

Бірліктер: операции умножения со сложением

Мына жағдайда орындалады

главные члены оценок; для семиточечной матрицы с ILU(0) $c$ — порядка двух-трёх десятков.

Прямой метод платит степенью размера задачи, итерационный — числом итераций, которое тоже растёт с размером, но медленнее.

Есептеу мысалыПрямое или итерационное решение для трёхмерной модели

Бастапқы деректер (шартты)

Сетка: 100 × 100 × 20 ячеек, по уравнению давления на ячейку: $n$ = 200 000
Итерационный решатель (допущение): 50 итераций сопряжённых градиентов с предобусловливанием, около 25 операций на неизвестное за итерацию
Компьютер (допущение): $10^9$ операций умножения со сложением в секунду на одно ядро
Число решений за расчёт (допущение): 4000: сотни шагов по несколько итераций

Шешім

Плотная матрица: $n^3/3 = 2{,}67\cdot10^{15}$ операций — 31 сут на одно решение.
Лента при естественной нумерации: $b = 100\cdot100 = 10\,000$ ; работа $nb^2 = 2\cdot10^{13}$ — 5,6 ч; память $nb = 2\cdot10^{9}$ чисел, по 8 байт — 16 ГБ.
Нумерация сначала по вертикали: $b = 20\cdot100 = 2000$ ; работа $8\cdot10^{11}$ — 13 мин на решение, память 3,2 ГБ.
Итерационный решатель: $50\cdot25\cdot200\,000 = 2{,}5\cdot10^{8}$ операций — 0,25 с.
За 4000 решений: лента (лучшая нумерация) — 37 сут, итерационный решатель — 17 мин.

Жауап: Ленточное исключение заняло бы 13 мин на одно решение и 37 сут на расчёт, итерационный решатель — 0,25 с и 17 мин.

Бұл нені білдіреді. Для трёхмерных моделей выбор очевиден, и он тем очевиднее, чем больше сетка: работа ленточного метода растёт как $n^{7/3}$ (ширина ленты — слой ячеек), итерационного — примерно как $n$ , умноженное на медленно растущее число итераций. Скорость компьютера и число операций здесь — допущения для оценки порядка величин. Данные примера условные.

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

Если решатель сходится медленно или не сходится

1
Читать журнал
на каких шагах растёт число линейных итераций и где сокращается шаг; после пуска и остановки скважин это нормально, на спокойных участках — сигнал.
2
Искать плохие ячейки
крошечный поровый объём рядом с крупными ячейками, проводимости, отличающиеся на много порядков, несоседние соединения с огромной проводимостью: они портят обусловленность; такие ячейки делают неактивными или объединяют с соседними.
3
Проверить скважины
большой индекс скважины в ячейке с малым объёмом даёт резкие изменения давления и плохую обусловленность.
4
Менять настройки осознанно
допуски и тип предобусловливателя меняют, сверяя баланс и результаты с расчётом при прежних настройках.

Терминдер мен белгілер

Линейный решатель: часть симулятора, решающая системы $\mathbf{A}\mathbf{x} = \mathbf{b}$ неявных шагов и итераций Ньютона.
Разреженная матрица: матрица, в строке которой ненулевых чисел столько, сколько связей у ячейки, плюс диагональ.
Ширина ленты $b$: наибольшая разность номеров связанных неизвестных; при естественной нумерации — число ячеек в строке (2D) или в слое (3D).
Заполнение: новые ненулевые элементы, которые появляются при исключении: исключённое неизвестное связывает своих соседей.
Спектральный радиус $\rho$: наибольший модуль собственного числа матрицы перехода итерационного метода; метод сходится при $\rho < 1$ .
Число обусловленности $\kappa$: отношение наибольшего собственного числа симметричной матрицы к наименьшему; во сколько раз относительная ошибка может превышать относительную невязку.
Предобусловливатель: легко обратимое подобие матрицы, уменьшающее число итераций: ILU, вложенная факторизация, CPR, многосеточные методы.
Невязка: $\mathbf{b} - \mathbf{A}\mathbf{x}$ ; для уравнения давления — небаланс ячейки в м³/сут.

Жиі қателіктер

Қате

«Чтобы решить систему, надо обратить матрицу».

Шын мәнінде

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

Қате

«Невязка мала — значит, давления точны».

Шын мәнінде

Относительная ошибка может превышать относительную невязку в $\kappa$ раз. В примере урока при невязке $10^{-2}$ изменение давления в ячейке 2 вышло почти вчетверо больше точного, а баланс шага нарушился на 1,9 % закачки.

Қате

«Если расчёт медленный, можно ослабить допуск решателя».

Шын мәнінде

Недорешённая система «создаёт» или «теряет» флюид, и ошибка баланса накапливается от шага к шагу; в полностью неявной схеме грубый допуск ещё и добавляет итераций Ньютона. Допуск меняют, только сверяя баланс и результаты.

Қате

«Простые итерации годятся для любой сетки — понадобится лишь больше итераций».

Шын мәнінде

Их число растёт как квадрат числа ячеек вдоль стороны: в модельной задаче 40 × 40 Якоби нужно 3729 итераций, а сопряжённым градиентам с ILU(0) — 30. Для больших моделей нужны методы Крылова с предобусловливанием.

Толық: тереңдетілген

Почему простые итерации медленны на больших сетках

Запишем матрицу как $\mathbf{A} = \mathbf{D} - \mathbf{L} - \mathbf{U}$ : диагональ, нижняя и верхняя части с обратным знаком. Метод Якоби: $\mathbf{D}\mathbf{x}^{(\nu+1)} = \mathbf{b} + (\mathbf{L} + \mathbf{U})\,\mathbf{x}^{(\nu)}$ . Вычтем его из точного равенства $\mathbf{D}\mathbf{x} = \mathbf{b} + (\mathbf{L} + \mathbf{U})\,\mathbf{x}$ :

Ошибка итераций и матрица перехода

\mathbf{e}^{(\nu+1)} = \mathbf{G}\,\mathbf{e}^{(\nu)},\qquad \mathbf{G}_J = \mathbf{D}^{-1}(\mathbf{L} + \mathbf{U}),\qquad \mathbf{G}_{GS} = (\mathbf{D} - \mathbf{L})^{-1}\,\mathbf{U}

мұндағы:

$\mathbf{e}^{(\nu)}$ ошибка приближения $\mathbf{x} - \mathbf{x}^{(\nu)}$ , в единицах неизвестных
$\mathbf{D},\ \mathbf{L},\ \mathbf{U}$ диагональ, нижняя и верхняя треугольные части матрицы (с обратным знаком), в единицах матрицы
$\mathbf{G}_J,\ \mathbf{G}_{GS}$ матрицы перехода методов Якоби и Зейделя, —

Бірліктер: безразмерные

Мына жағдайда орындалады

линейная система с постоянной матрицей; диагональ без нулей.

После $\nu$ итераций ошибка равна $\mathbf{G}^{\nu}\mathbf{e}^{(0)}$ : её составляющие вдоль собственных векторов $\mathbf{G}$ умножаются на собственные числа. Медленнее всех затухает составляющая с наибольшим модулем собственного числа — спектральным радиусом.

Для одномерной полосы из $N$ одинаковых ячеек с заданным давлением за концами, в простейшей записи (на диагонали 2, рядом −1), собственные векторы матрицы Якоби — синусоиды с $m = 1, \ldots, N$ полуволнами, а собственные числа — $\cos[m\pi/(N + 1)]$ . Самая гладкая синусоида затухает медленнее всех: $\rho_J = \cos[\pi/(N + 1)] \approx 1 - \pi^2/[2(N + 1)^2]$ , и число итераций $\nu \approx 2(N + 1)^2\ln(1/\varepsilon)/\pi^2$ — пропорционально квадрату числа ячеек. Изрезанные составляющие, меняющие знак от ячейки к ячейке, метод Зейделя убирает за несколько итераций; на этом построены многосеточные методы: простые итерации сглаживают ошибку, а её гладкую часть исправляют на грубой сетке, где она уже не гладкая.

Как устроен метод сопряжённых градиентов

Алгоритм метода сопряжённых градиентов

\begin{aligned} &\mathbf{r}^{(0)} = \mathbf{b} - \mathbf{A}\mathbf{x}^{(0)},\quad \mathbf{p}^{(0)} = \mathbf{r}^{(0)};\\ &\alpha_\nu = \frac{\mathbf{r}^{(\nu)\mathsf T}\mathbf{r}^{(\nu)}}{\mathbf{p}^{(\nu)\mathsf T}\mathbf{A}\,\mathbf{p}^{(\nu)}},\quad \mathbf{x}^{(\nu+1)} = \mathbf{x}^{(\nu)} + \alpha_\nu\,\mathbf{p}^{(\nu)},\quad \mathbf{r}^{(\nu+1)} = \mathbf{r}^{(\nu)} - \alpha_\nu\,\mathbf{A}\,\mathbf{p}^{(\nu)};\\ &\beta_\nu = \frac{\mathbf{r}^{(\nu+1)\mathsf T}\mathbf{r}^{(\nu+1)}}{\mathbf{r}^{(\nu)\mathsf T}\mathbf{r}^{(\nu)}},\quad \mathbf{p}^{(\nu+1)} = \mathbf{r}^{(\nu+1)} + \beta_\nu\,\mathbf{p}^{(\nu)}\end{aligned}

мұндағы:

$\mathbf{r}^{(\nu)}$ невязка на итерации $\nu$ , в единицах правой части
$\mathbf{p}^{(\nu)}$ направление поиска, в единицах неизвестных
$\alpha_\nu,\ \beta_\nu$ длина шага вдоль направления и коэффициент нового направления, —

Бірліктер: в единицах системы

Мына жағдайда орындалады

симметричная положительно определённая матрица; с предобусловливателем в направлениях и скалярных произведениях вместо $\mathbf{r}$ используют $\mathbf{z} = \mathbf{M}^{-1}\mathbf{r}$ .

За итерацию — одно умножение матрицы на вектор и два скалярных произведения. Направления сопряжены: $\mathbf{p}^{(\mu)\mathsf T}\mathbf{A}\,\mathbf{p}^{(\nu)} = 0$ при $\mu \ne \nu$ , поэтому шаг вдоль нового направления не портит достигнутого вдоль прежних. В точной арифметике после $k$ итераций ошибка наименьшая из возможных в подпространстве Крылова размерности $k$ ; машинные округления постепенно нарушают сопряжённость, и на практике метод работает как итерационный.

Системы полностью неявной схемы и двухступенчатое предобусловливание

В полностью неявной схеме у ячейки два-три неизвестных — давление и насыщенности, и матрица Ньютона состоит из блоков 2 × 2 или 3 × 3 (якобиан). Части задачи ведут себя по-разному. Давление ведёт себя почти «эллиптически»: на неявном шаге изменение дебита скважины сказывается во всей модели, и ошибка давления гладкая и общая. Насыщенности «гиперболические»: они переносятся потоком и меняются в основном у фронтов и скважин, и их ошибка местная.

Двухступенчатое предобусловливание CPR использует это различие. На первой ступени из уравнений фаз ячейки составляют одно уравнение для давления — примерно так, как суммируют балансы фаз в уравнении давления (уравнение давления), — и решают его сильным методом, обычно алгебраическим многосеточным. На второй ступени дешёвое неполное LU-разложение всей системы убирает оставшуюся местную ошибку. Каждая ступень делает то, что умеет лучше, и число итераций внешнего метода (как правило, GMRES) остаётся небольшим даже для моделей в миллионы ячеек.

Басқа тақырыптармен байланыс

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

Өзіңізді тексеріңіз

1
Сколько ненулевых чисел в строке матрицы давления для внутренней ячейки трёхмерной декартовой сетки без несоседних соединений и скважин?
1. А.7
2. Б.3
3. В.27
4. Г.5
Жауабын көрсету
Дұрыс жауап: А. 7
Диагональ и шесть соседей по граням — справа, слева, спереди, сзади, сверху и снизу. 27 чисел было бы, если связывать ячейку со всеми соседями по углам и рёбрам, как в некоторых многоточечных схемах.
2
Модель — 200 × 100 × 30 ячеек. Чему равна ширина ленты при естественной нумерации (быстрее всего — вдоль длинной оси, затем вдоль второй, затем по слоям) и при нумерации сначала по вертикали, затем вдоль короткой горизонтальной оси? Во сколько раз различается работа ленточного исключения?
Жауабын көрсету
Естественная нумерация: $b = 200\cdot100 = 20\,000$ . Сначала по вертикали, затем вдоль короткой оси: $b = 30\cdot100 = 3000$ . Работа пропорциональна $b^2$ : $(20\,000/3000)^2 \approx 44$ . Но даже 3000 — слишком широкая лента для шести сотен тысяч неизвестных: $nb^2 \approx 5{,}4\cdot10^{12}$ операций на одно решение.
3
Шаг по времени увеличили в 10 раз. Как изменится число обусловленности системы давления закрытого пласта со скважинами с заданным дебитом, если шаг уже крупный?
1. А.Не изменится
2. Б.Уменьшится примерно в 10 раз
3. В.Увеличится в 100 раз
4. Г.Увеличится примерно в 10 раз
Жауабын көрсету
Дұрыс жауап: Г. Увеличится примерно в 10 раз
Наименьшее собственное число равно $S/\Delta t$ — оно уменьшится в 10 раз, а наибольшее почти не изменится: его задают связи. В примере урока при переходе от 10 к 100 сут число обусловленности выросло с 14,6 до 137, а итераций Якоби понадобилось в семь раз больше.
4
Методу Якоби на сетке 30 × 30 понадобилось 600 итераций, а сопряжённым градиентам — 50. Сколько примерно понадобится каждому методу на сетке 60 × 60 того же пласта?
Жауабын көрсету
Число итераций Якоби растёт как квадрат числа ячеек вдоль стороны: около $600\cdot4 = 2400$ . У сопряжённых градиентов — как корень из числа обусловленности, то есть пропорционально числу ячеек вдоль стороны: около 100. При этом каждая итерация на сетке 60 × 60 вчетверо дороже, так что общая работа растёт в 16 и в 8 раз.
5
Относительная невязка решения — $10^{-3}$ , число обусловленности системы — около $10^4$ . Что можно сказать об относительной ошибке давлений?
1. А.Она не больше $10^{-7}$
2. Б.Она может достигать $10$ , то есть быть больше самого решения
3. В.Она не больше $10^{-3}$
4. Г.Она равна нулю, раз невязка мала
Жауабын көрсету
Дұрыс жауап: Б. Она может достигать $10$ , то есть быть больше самого решения
Оценка сверху — $\kappa\cdot10^{-3} = 10$ : при плохой обусловленности малая невязка не гарантирует даже верного порядка величин. На практике ошибка обычно меньше оценки, но поэтому в плохо обусловленных задачах допуск делают строже или улучшают предобусловливатель.
6
Сумма невязок уравнений давления по всем ячейкам после решения — +15 м³/сут при шаге 20 сут. Что это значит для материального баланса?
Жауабын көрсету
Сумма невязок — небаланс модели: перетоки между ячейками при сложении уничтожаются, и остаются накопление и скважины. За шаг уравнения давления «потеряли» или «создали» $15\cdot20 = 300$ м³ флюида — в зависимости от знака, принятого в записи невязки. Это и есть погрешность материального баланса, которую печатает симулятор; её сравнивают с отборами и поровым объёмом и при заметной величине ужесточают допуск.

Жаңартылған: 19 октября 2026 г.