\chapter{The Newton-Raphson method}

For the isothermal immiscible multi-phase model, the following mass
conservation equation needs to be solved:
\begin{align}
\underbrace{
  \frac{\partial \phi \varrho_\alpha S_\alpha}{\partial t}
 -
 \text{div} \left\{
 \varrho_\alpha \frac{k_{r\alpha}}{\mu_\alpha} \mbox{\textbf{K}} \left(\text{grad}\, p_\alpha - \varrho_{\alpha} \mbox{\textbf{g}} \right)
 \right\} - q_\alpha} _
{\textbf{f}(\textbf{u})}
= 0 \; .
\end{align}

Because of the nonlinear dependencies even solving this comparatively
simple equation is a quite challenging task. However, for
finding roots of non-linear systems equations the
\textsc{Newton}-\textsc{Raphson} method can be used.

When using a fully-implicit numerical model, each time step essentially
consists of the application of the \textsc{Newton} algorithm to solve
the nonlinear system.

The idea of this algorithm is to linearize the non-linear system of
equation at a given solution, and then solving the resulting linear
system of equations. The hope is, that the solution of this linear
system is closer to the root of the non-linear system of
equations. This process is repeated until either convergence is
reached (a pre-determined accuracy is reached), or divergence of the
algorithm is detected (either by trespassing the maximum number of
iterations or by failure to linearize). This method can be formalized
as follows:

\begin{subequations}
\begin{align}
\label{NewtonGen}
\textbf{u}^{r+1} &= \textbf{u}^r + \Delta \textbf{u}^r \\
\Delta \textbf{u}^r & = - \left\{\text{grad}\,\textbf{f} (\textbf{u}^r) \right\}^{-1} \textbf{f}(\textbf{u}^r) \\
\end{align}
\end{subequations}

\noindent with
\begin{itemize}
\item $\textbf{u}$: Vector of unknowns
\item $\textbf{f}(\textbf{u}^r)$: Residual (Function of the vector of unknowns which ought to be set to $0$)
\item $\phantom{a}^r$: last iteration, $\phantom{a}^{r+1}$: current iteration,
\item $\text{grad}\,\phantom{a}$: \textsc{Jacobian} matrix of
  $\textbf{f}$, i.e. matrix of the derivatives of \textbf{f} regarding
  all components of $\textbf{u}$
\end{itemize}

The value of $\textbf{u}$ for which $\textbf{f}$ becomes zero is
searched for. Bringing \eqref{NewtonGen} into the form used the linear
solvers
\begin{equation}
\label{GenSysEq}
 \textbf{A}\textbf{x} - \textbf{b} = 0 
\end{equation}
leads to
\begin{itemize}
\item $\textbf{A} = \text{grad}\,\textbf{f} (\textbf{u}^r)$
\item $\textbf{x} = \textbf{u}^{r} - \textbf{u}^{r+1}$
\item $\textbf{b} = \textbf{f}(\textbf{u}^{r})$
\end{itemize}

Once $\textbf{u}^{r} - \textbf{u}^{r+1}$ has been calculated, \eWoms
updates the current solution in \texttt{NewtonController::update()}
and starts the next iteration if the scheme has not yet converged.

%%% Local Variables: 
%%% mode: latex
%%% TeX-master: "ewoms-handbook"
%%% End: