Basic numerics

In this section we discuss some basic numerical concepts for the solution of partial differential equations used in DuMux.

Newton's method

Nonlinear systems of equations can be solved with Newton's method. Newton's method is quadratically convergent, if the initial solution is close enough to the root of the residual equation. It does however not guarantee global convergence.

The discrete differential equations are formulated in residual form. All terms are on the left hand side and are summed up. The terms contain values for the primary variables which are part of the solution vector $\mathbf{\text{u}}$. The sum of the terms is called residual $\mathbf{\text{r}}(\mathbf{\text{u}})$ which is a function of the solution. For example:

$$\begin{array}{r}\underset{=:\mathbf{\text{r}}(\mathbf{\text{u}})}{\underbrace{\varphi \frac{\partial {\varrho }_{\alpha }{S}_{\alpha }}{\partial t}-\text{div}\left({\varrho }_{\alpha }\frac{k_{r\alpha }}{{\mu }_{\alpha }}\mathit{K}(\mathrm{\nabla }{p}_{\alpha }-{\varrho }_{\alpha }\mathit{g})\right)-q_{\alpha }}}=0\end{array}$$

We are looking for the unknown solution $\mathbf{\text{u}}$ and use Newton's method to create a series of ${\mathbf{\text{u}}}^i$ that we hope converges to $\mathbf{\text{u}}$. We start with an initial guess ${\mathbf{\text{u}}}^0$ and calculate the initial residual $\mathbf{\text{r}}({\mathbf{\text{u}}}^0)$. Next, we calculate the derivative of the residual with respect to the solution. This is the Jacobian matrix

$$\begin{array}{r}\frac{\text{d}}{\text{d}\mathbf{\text{u}}}\mathbf{\text{r}}\left({\mathbf{\text{u}}}^i\right)=J_{\mathbf{\text{r}}\left({\mathbf{\text{u}}}^i\right)}={\left(\frac{\text{d}}{\text{d}{\mathbf{\text{u}}}_m^i}\mathbf{\text{r}}{\left({\mathbf{\text{u}}}^i\right)}_n\right)}_{m,n}\end{array}$$

with $i$ denoting the Newton iteration step. Each column is the residual derived with respect to the $m$th entry of ${\mathbf{\text{u}}}^i$.

The Jacobian indicates the direction where the residual increases. By solving the linear system

$$\begin{array}{r}{J}_{\mathbf{\text{r}}({\mathbf{\text{u}}}^i)}\cdot {\mathbf{\text{x}}}^i=\mathbf{\text{r}}({\mathbf{\text{u}}}^i)\end{array}$$

we calculate the direction of maximum growth ${\mathbf{\text{x}}}^i$. We subtract it from our current solution to get a new, better solution ${\mathbf{\text{u}}}^{i+1}={\mathbf{\text{u}}}^i-{\mathbf{\text{x}}}^i$.

We repeat the calculation of the Jacobian $J_{\mathbf{\text{r}}({\mathbf{\text{u}}}^i)}$ and the direction of maximum growth ${\mathbf{\text{x}}}^i$ until our approximated solution becomes good enough. See Dumux::NewtonSolver for the implementation of the Newton method based solver in DuMux which features various convergence criteria and a simple line search algorithm to improve the update.

Time discretization

In the following, we describe two basic time stepping schemes which are first-order accurate. Our systems of partial differential equations are discretized in space and in time.

Let us consider the general case of a balance equation of the following form

$$\frac{\partial m(u)}{\partial t}+\mathrm{\nabla }\cdot \mathit{f}(u,\mathrm{\nabla }u)+q(u)=0,$$

seeking an unknown quantity $u$ in terms of storage $m$, flux $\mathit{f}$ and source $q$. All available Dumux models can be written mathematically in this form with possibly vector-valued quantities $u$, $m$, $q$ and a tensor-valued flux $\mathit{f}$. For the sake of simplicity, we here assume scalar quantities $u$, $m$, $q$ and a vector-valued flux $\mathit{f}$.

To discretize the continuous form of the balance equations, we need to choose an approximation for the temporal derivative $\partial m(u)/\partial t$. While many elaborate methods for this approximation exist, we focus on the simplest one of a first order difference quotient

$$\frac{\partial m(u_{k/k+1})}{\partial t}\approx \frac{m(u_{k+1})-m(u_k)}{\mathrm{\Delta }{t}_{k+1}}$$

for approximating the solution $u$ at time $t_k$ (forward) or $t_{k+1}$ (backward). The question of whether to choose the forward or the backward quotient leads to the explicit and implicit Euler method, respectively. Using the explicit Euler method leads to

$$\frac{m(u_{k+1})-m(u_k)}{\mathrm{\Delta }{t}_{k+1}}+\mathrm{\nabla }\cdot \mathit{f}(u_k,\mathrm{\nabla }{u}_k)+q(u_k)=0,$$

whereas the implicit Euler method is described as

$$\frac{m(u_{k+1})-m(u_k)}{\mathrm{\Delta }{t}_{k+1}}+\mathrm{\nabla }\cdot \mathit{f}(u_{k+1},\mathrm{\nabla }{u}_{k+1})+q(u_{k+1})=0.$$

Once the solution $u_k$ at time $t_k$ is known, it is straightforward to determine $m(u_{k+1})$, while attempting to do the same for the equation discreized with the implicit Euler involves the solution of a system of equations. The explicit method is stable only if the time step size $\mathrm{\Delta }{t}_{k+1}$ is below a certain limit that depends on the specific balance equation, whereas the implicit method is unconditionally stable.

---

🛠 Edit above doc on GitLab