Методы решения нелинейных адвекционно-диффузионных систем за пределами Ньютона-Рафсона?

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

\frac{\partial с_{1}}{\partial T} знак равно 0 знак равно F_{1} + Q_{1} (с_{1}, с_{2}) \frac{\partial с_{2}}{\partial T} знак равно 0 знак равно F_{2} + Q_{2} (с_{1}, с_{2})

$\frac{\partial c_1}{\partial t} = 0 = \mathcal{F}_1 + \mathcal{Q}_1(c_1,c_2) \\ \frac{\partial c_2}{\partial t} = 0 = \mathcal{F}_2 + \mathcal{Q}_2(c_1,c_2)$

куда $\mathcal{F}_i$ диффузный и адвективный поток для видов $i$ , а также $\mathcal{Q}_i$ является исходным термином для вида $i$ ,

Я смог написать решатель для моей проблемы, используя метод Ньютона-Рафсона, и полностью связал две области, используя матрицу блочной массы, то есть:

F_{с о U п L е d} знак равно [\begin{array}{cc} A_{1} & 0 \\ 0 & A_{2} \end{array}] \underset{{Икс}_{я}}{\underset{⏟}{[\begin{matrix} с_{1, я} \\ с_{2, я} \end{matrix}]}} - [\begin{matrix} б_{1} (с_{1, я}, с_{2, я}) \\ б_{2} (с_{1, я}, с_{2, я}) \end{matrix}]

$F_{coupled} = \left[\begin{array}{c c} A_1 & 0 \\ 0 & A_2 \\ \end{array}\right]\underbrace{ \left[\begin{array}{c} c_{1,i} \\ c_{2,i} \\ \end{array}\right] }_{x_i} - \left[\begin{array}{c} b_1(c_{1,i}, c_{2,i}) \\ b_2(c_{1,i}, c_{2,i}) \\ \end{array}\right]$

Семестр $F_{coupled}$ используется для определения матрицы Якоби и обновления как $c_1$ а также $c_2$ :

J ({Икс}_{я}) [{Икс}_{я + 1} - {Икс}_{я}] знак равно - F_{с о U п L е d}

$\mathcal{J}(x_i) \left[x_{i+1} - x_i \right] = -\mathcal{F}_{coupled}$

или

{Икс}_{я + 1} знак равно {Икс}_{я} - {(J ({Икс}_{я}))}^{- 1} F_{с о U п L е d}

$x_{i+1} = x_i - \left(\mathcal{J}(x_i)\right)^{-1} \mathcal{F}_{coupled}$

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

Проблема в том, что я собираюсь перейти к более крупной системе, в которой оба домена находятся в 2D / 2.5D, и вычисление матрицы Якоби приведет к быстрому истощению моих доступных компьютерных ресурсов. Я строю эту модель для использования в настройках оптимизации позже, поэтому я не могу сидеть за рулем на каждой итерации, настраивая коэффициент демпфирования и т. Д.

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

Могут ли я упустить какие-либо другие хитрые алгоритмы, которые могут решить систему нелинейных уравнений, не прибегая к перерасчету якобиана как обидного?

nonlinear-equations newton-method advection-diffusion coupling cbcoutinho
источник

Рассматривали ли вы итеративные решатели, такие как AMG - алгебраические многосеточные методы. Возможно, вам придется придумать хорошие предварительные кондиционеры, основанные на физике.

NameRakes

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

Билл Барт

Нет, я не рассматривал AMG, я думал, что они предназначены только для симметричных систем и не могут быть использованы в задачах с преобладанием конвекции. Я снова посмотрю в литературе для AMG.

cbcoutinho

Параллельные вычисления сложны, потому что этот проект разрабатывается как отдельное приложение для коллег, у которых нет доступа к таким ресурсам. Я подумывал о том, чтобы превратить mpi в проект для себя, но это увеличило бы барьер входа для других, что и было

главной целью

Почему вычисление якобиана так проблематично? Если вы используете конечные различия / объемы / элементы, он должен иметь разреженную часть, которая всегда одинакова, и диагональную часть, которая изменяется, но вычисляется тривиально.

Дэвид Кетчесон

Ответы:

Я предполагаю, что ограничение в 2D и 3D хранит якобиан.

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

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

D F (U^{N}) δ U^{N} знак равно - F (U^{N})

$DF(u^n)\, \delta u^n = -F(u^n)$ (где

D F

$DF$ Якобиан) можно решить с помощью решателя подпространств Крылова, такого как GMRES и BiCGStab, используя тот факт, что

D F (U^{N}) δ U \approx \frac{F (U^{N} + ε \frac{δ U}{| | δ U | |}) - F (U^{N})}{ε},

$DF(u^n)\, \delta u \approx \frac{F(u^n+\epsilon\frac{\delta u}{\Vert\delta u\Vert}) - F(u^n)}{\epsilon}.$ Это потому, что GMRES и BiCGStab не требуют матрицы LHS

A

$A$ им нужно только вычислить свой продукт

A x

$Ax$ заданный вектор

x

$x$ ,

Теперь с надлежащим значением $\epsilon$ (обычно о $10^{-7}$ для операций с плавающей запятой двойной точности) вы можете выполнить цикл Ньютона, не вычисляя и не сохраняя якобиан. Я точно знаю, что эта техника используется для решения некоторых нетривиальных случаев в вычислительной гидродинамике. Обратите внимание, однако, что количество оценок функции $F$ будет больше, чем в методе хранения матрицы, вместо того, чтобы требовать матрично-векторного произведения.

Еще одна вещь, на которую следует обратить внимание: если ваша система такова, что необходим мощный прекондиционер (т.е. Jacobi или block-Jacobi не хватит), вы можете попробовать использовать вышеупомянутый метод в качестве сглаживающего в многосеточной схеме. Если вы хотите попробовать прекондиционер Якоби с точечной или блок-схемой, вы можете вычислить и сохранить только диагональные элементы или диагональные блоки якобиана, что немного. Я также хотел бы упомянуть, что предварительный кондиционер Гаусса-Зейделя или SSOR можно реализовать без явного сохранения якобиана. В данной статье описывается реализация безматричной GMRES, предварительно обусловленной безматричной симметричной Гаусса-Зейделя, в контексте вычислительной гидродинамики.

Адитья Каши
источник

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

Если вам просто нужна быстрая числовая схема для решения эволюции времени, взгляните на IMEX (неявно-явные) схемы; см., например, эту статью Ашера, Руута, Неявно-Явные Спитери Методы Рунге-Кутты для зависящих от времени дифференциальных уравнений в частных производных .

Вы можете даже попытаться использовать явную схему интеграции времени высокого порядка с контролем размера шага (как в Matlab ODE45). Тем не менее, вы можете столкнуться с проблемами из-за жесткости вашей системы, которая исходит от диффузионной части. К счастью, диффузионная часть является линейной, так что ее можно трактовать неявно (такова идея схем IMEX).

январь
источник

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

Во-первых, глядя на вашу запись $F_{coupled}$ Я не вижу связанную форму, я полагаю, что $b_1$ а также $b_2$ оба должны быть зависимы от $c_{1,i}$ а также $c_{2,i}$ , Более того, если $A_1$ а также $A_2$ являются матричными представлениями приближений $\mathcal{F}_1$ а также $\mathcal{F}_2$ тогда они не должны зависеть только от $c_i$ , но также и о соседских ценностях, но это может быть только неправильное понимание вашей записи.

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

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

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

Петр Фролкович
источник

Привет @ Питер, ты прав насчет связи, я отредактировал главное уравнение, чтобы показать и

b_{1}

$b_1$ а также

b_{2}

$b_2$ обе функции

c_{1}

$c_1$ а также

c_{2}

$c_2$ , Матрицы

A_{1}

$A_1$ а также

A_{2}

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

F_{1}

$F_1$ а также

F_{2}

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

cbcoutinho