From f24df1a6e1512df5395eca2899595f0ccff14f26 Mon Sep 17 00:00:00 2001
From: Timo Koch <timokoch@uio.no>
Date: Thu, 29 Jun 2023 00:31:06 +0200
Subject: [PATCH] [doxygen] Add numerics section from handbook
---
doc/doxygen/DoxygenDumuxLayout.xml | 1 +
doc/doxygen/pages/numerics.md | 113 +++++++++++++++++++
doc/handbook/0_dumux-handbook.tex | 6 -
doc/handbook/5_assemblinglinearsystem.tex | 123 ---------------------
doc/handbook/6_temporaldiscretizations.tex | 66 -----------
5 files changed, 114 insertions(+), 195 deletions(-)
create mode 100644 doc/doxygen/pages/numerics.md
delete mode 100644 doc/handbook/5_assemblinglinearsystem.tex
delete mode 100644 doc/handbook/6_temporaldiscretizations.tex
diff --git a/doc/doxygen/DoxygenDumuxLayout.xml b/doc/doxygen/DoxygenDumuxLayout.xml
index b8dcb4f257..40aa7c4931 100644
--- a/doc/doxygen/DoxygenDumuxLayout.xml
+++ b/doc/doxygen/DoxygenDumuxLayout.xml
@@ -25,6 +25,7 @@ SPDX-License-Identifier: GPL-3.0-or-later
</tab>
<tab type="modules" visible="yes" title="Module documentation" intro=""/>
<tab type="user" url="@ref flow-and-transport-in-porous-media" visible="yes" title="Flow and transport in porous media" intro=""/>
+ <tab type="user" url="@ref basic-numerics" visible="yes" title="Basic numerics" intro=""/>
<tab type="user" url="@ref python-bindings" visible="yes" title="Python bindings" intro=""/>
<tab type="usergroup" title="Developer documentation">
<tab type="user" url="@ref dumux-and-dune" visible="yes" title="Dumux and Dune" intro=""/>
diff --git a/doc/doxygen/pages/numerics.md b/doc/doxygen/pages/numerics.md
new file mode 100644
index 0000000000..809e8aa47f
--- /dev/null
+++ b/doc/doxygen/pages/numerics.md
@@ -0,0 +1,113 @@
+# 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 $\textbf{u}$. The sum of the terms
+is called residual $\textbf{r}(\textbf{u})$ which is a function of the solution. For
+example:
+
+\begin{align*}
+\underbrace{
+ \phi \frac{\partial \varrho_\alpha S_\alpha}{\partial t}
+ -
+ \text{div} \left(
+ \varrho_\alpha \frac{k_{r\alpha}}{\mu_\alpha} \boldsymbol{K}
+ \left(\nabla p_\alpha - \varrho_{\alpha} \boldsymbol{g} \right)
+ \right) - q_\alpha} _
+{=: \, \textbf{r}(\textbf{u})}
+= 0
+\end{align*}
+
+We are looking for the unknown solution $\textbf{u}$ and use Newton's method to
+create a series of $\textbf{u}^i$ that we hope converges to $\textbf{u}$.
+We start with an initial guess $\textbf{u}^0$ and
+calculate the initial residual $\textbf{r}(\textbf{u}^0)$. Next,
+we calculate the derivative of the residual with respect to the solution.
+This is the Jacobian matrix
+
+\begin{align*}
+ \frac{\text{d}}{\text{d}\textbf{u}}\textbf{r} \left(\textbf{u}^i\right)
+ = J_{\textbf{r} \left(\textbf{u}^i\right)}
+ = \left(\frac{\text{d}}{\text{d}\textbf{u}^i_m}\textbf{r} \left(\textbf{u}^i\right)_n\right)_{m,n}
+\end{align*}
+
+with $i$ denoting the Newton iteration step.
+Each column is the residual derived with respect to the $m$th entry of $\textbf{u}^i$.
+
+The Jacobian indicates the direction where the residual increases. By solving the
+linear system
+
+\begin{align*}
+ J_{\textbf{r}(\textbf{u}^i)} \cdot \textbf{x}^i = \textbf{r}(\textbf{u}^i)
+\end{align*}
+
+we calculate the direction of maximum growth $\textbf{x}^i$. We subtract it from
+our current solution to get a new, better solution
+$\textbf{u}^{i+1} = \textbf{u}^i - \textbf{x}^i$.
+
+We repeat the calculation of the Jacobian $J_{\textbf{r}(\textbf{u}^i)}$ and the
+direction of maximum growth $\textbf{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
+
+\begin{equation}
+\frac{\partial m(u)}{\partial t} + \nabla\cdot\boldsymbol{f}(u, \nabla u) + q(u) = 0,
+\end{equation}
+
+seeking an unknown quantity $u$ in terms of storage $m$, flux $\boldsymbol{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 $\boldsymbol{f}$.
+For the sake of simplicity, we here assume scalar quantities $u$, $m$, $q$ and a vector-valued
+flux $\boldsymbol{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
+
+\begin{equation}
+\frac{\partial m(u_{k/k+1})}{\partial t}
+\approx \frac{m(u_{k+1}) - m(u_k)}{\Delta t_{k+1}}
+\end{equation}
+
+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
+
+\begin{equation}
+\frac{m(u_{k+1}) - m(u_k)}{\Delta t_{k+1}} + \nabla\cdot\boldsymbol{f}(u_k, \nabla u_k) + q(u_k) = 0,
+\end{equation}
+
+whereas the implicit Euler method is described as
+
+\begin{equation}
+\frac{m(u_{k+1}) - m(u_k)}{\Delta t_{k+1}}
++ \nabla\cdot\boldsymbol{f}(u_{k+1}, \nabla u_{k+1}) + q(u_{k+1}) = 0.
+\end{equation}
+
+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 $\Delta t_{k+1}$ is below a certain limit that depends
+on the specific balance equation, whereas the implicit method
+is unconditionally stable.
diff --git a/doc/handbook/0_dumux-handbook.tex b/doc/handbook/0_dumux-handbook.tex
index be940e0112..a90966bcc6 100644
--- a/doc/handbook/0_dumux-handbook.tex
+++ b/doc/handbook/0_dumux-handbook.tex
@@ -120,12 +120,6 @@ The handbook is in the process of being dissolved.
\input{5_structure}
\input{5_newfoldersetup}
\input{5_scripts}
-\input{5_assemblinglinearsystem}
-
-\chapter{Advanced \Dumux\ -- Detailed Instructions}
-This chapter contains detailed information for those who are interested
-in deeper modifications of underlying \Dumux models, classes, functions, etc.
-\input{6_temporaldiscretizations}
\bibliographystyle{plainnat}
\bibliography{dumux-handbook}
diff --git a/doc/handbook/5_assemblinglinearsystem.tex b/doc/handbook/5_assemblinglinearsystem.tex
deleted file mode 100644
index affcf49ba3..0000000000
--- a/doc/handbook/5_assemblinglinearsystem.tex
+++ /dev/null
@@ -1,123 +0,0 @@
-% SPDX-FileCopyrightInfo: Copyright © DuMux Project contributors, see AUTHORS.md in root folder
-% SPDX-License-Identifier: CC-BY-4.0
-
-\section{Assembling the linear system}
-\label{sc_linearsystem}
-The physical system is implemented as the mathematical differential equation in
-local operators. \Dumux generates the linear system automatically. Read on, to
-learn what is done internally.
-
-\subsection{Newton's method}
-The differential equations are implemented in the 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 $\textbf{u}$. The sum of the terms
-is called residual $\textbf{r}(\textbf{u})$ which is a function of the solution. For
-example:
-\begin{align*}
-\underbrace{
- \phi \frac{\partial \varrho_\alpha S_\alpha}{\partial t}
- -
- \text{div} \left(
- \varrho_\alpha \frac{k_{r\alpha}}{\mu_\alpha} \mathbf{K}
- \left(\grad\, p_\alpha - \varrho_{\alpha} \mathbf{g} \right)
- \right) - q_\alpha} _
-{=: \, \textbf{r}(\textbf{u})}
-= 0
-\end{align*}
-
-We don't know the solution $\textbf{u}$, so we use the iterative Newton's method to
-obtain a good estimate of $\textbf{u}$. We start with an initial guess $\textbf{u}^0$ and
-calculate it's residual $\textbf{r}(\textbf{u}^0)$. To minimize the error, we calculate
-the derivative of the residual with respect to the solution. This is the Jacobian
-matrix
-\begin{align*}
- \frac{\text{d}}{\text{d}\textbf{u}}\textbf{r} \left(\textbf{u}^i\right)
- = J_{\textbf{r} \left(\textbf{u}^i\right)}
- = \left(\frac{\text{d}}{\text{d}\textbf{u}^i_m}\textbf{r} \left(\textbf{u}^i\right)_n\right)_{m,n}
-\end{align*}
-with $i$ denoting the Newton iteration step.
-Each column is the residual derived with respect to the $m$th entry of $\textbf{u}^i$.
-
-The Jacobian indicates the direction where the residual increases. By solving the
-linear system
-\begin{align*}
- J_{\textbf{r}(\textbf{u}^i)} \cdot \textbf{x}^i = \textbf{r}(\textbf{u}^i)
-\end{align*}
-we calculate the direction of maximum growth $\textbf{x}^i$. We subtract it from
-our current solution to get a new, better solution
-$\textbf{u}^{i+1} = \textbf{u}^i - \textbf{x}^i$.
-
-We repeat the calculation of the Jacobian $J_{\textbf{r}(\textbf{u}^i)}$ and the
-direction of maximum growth $\textbf{x}^i$ until our approximated solution becomes good enough.
-
-\subsection{Structure of matrix and vectors}
-To understand the meaning of an entry in the matrix or the vector of the linear system, we have
-to define their structure. Both have a block structure. Each block contains the degrees of
-freedom (also called variables or unknowns) for a control volume. The equation index is used
-to order the degrees of freedom. For each control volume we have one block. The mapper is
-used to order the blocks.
-
-\begin{figure}[htbp]
-\begin{center}
-\begin{tikzpicture}[fill=dumuxBlue]
- %% blocking structure
- % matrix
- \node at (0.3,4.2){\footnotesize 1. CV};
- \node at (1.7,4.2){\footnotesize 2. CV};
- \node at (3.5,4.2){\footnotesize $n$. CV};
-
- \draw (0,0) rectangle (4,4);
-
- \fill (0.1,3.1) rectangle (0.9,3.9);
- \fill (1.1,3.1) rectangle (1.9,3.9);
- \node at (2.5,3.5) {$\dots$};
- \fill (3.1,3.1) rectangle (3.9,3.9);
- \node at (4,3.5) [right]{\footnotesize 1. CV};
-
- \fill (0.1,2.1) rectangle (0.9,2.9);
- \fill (1.1,2.1) rectangle (1.9,2.9);
- \node at (2.5,2.5) {$\dots$};
- \fill (3.1,2.1) rectangle (3.9,2.9);
- \node at (4,2.5) [right]{\footnotesize 2. CV};
-
- \node at (0.5,1.5) {$\vdots$};
- \node at (1.5,1.5) {$\vdots$};
- \node at (2.5,1.5) {$\ddots$};
- \node at (3.5,1.5) {$\vdots$};
-
- \fill (0.1,0.1) rectangle (0.9,0.9);
- \fill (1.1,0.1) rectangle (1.9,0.9);
- \node at (2.5,0.5) {$\dots$};
- \fill (3.1,0.1) rectangle (3.9,0.9);
- \node at (4,0.5) [right]{\footnotesize $n$. CV};
-
- % vector
- \draw (5.5,0) rectangle (5.9,4);
- \fill (5.6,3.1) rectangle (5.8,3.9);
- \fill (5.6,2.1) rectangle (5.8,2.9);
- \node at (5.7,1.5) {$\vdots$};
- \fill (5.6,0.1) rectangle (5.8,0.9);
-
- %% intra-block structure
- \fill (8.1,2.1) rectangle (8.9,2.9);
- \draw (9,2.8) -- (9.6,3.4);
- \draw (9,2.6) -- (9.6,2.8);
- \draw (9,2.2) -- (9.3,1.6);
-
- \node at (10,4) {${eqIdx}$};
- \node at (10,3.4) {$0$};
- \node at (10,2.8) {$1$};
- \node at (10,2.2) {$\vdots$};
- \node at (10,1.6) {$m-1$};
-
- \fill (11.1,2.1) rectangle (11.3,2.9);
- \draw (11,2.8) -- (10.4,3.4);
- \draw (11,2.6) -- (10.4,2.8);
- \draw (11,2.2) -- (10.7,1.6);
-\end{tikzpicture}
-\end{center}
-\caption{Structure of matrix and vector, left: block structure, right: within block}
-\end{figure}
-
-Accessing entries follows this structure. You can access the pressure value in the third ($n=3$) sub-control volume in
-a solution vector \lstinline{sol} with \lstinline{sol[n-1][pressureIdx]=sol[2][pressureIdx]}.
diff --git a/doc/handbook/6_temporaldiscretizations.tex b/doc/handbook/6_temporaldiscretizations.tex
deleted file mode 100644
index 5608500442..0000000000
--- a/doc/handbook/6_temporaldiscretizations.tex
+++ /dev/null
@@ -1,66 +0,0 @@
-% SPDX-FileCopyrightInfo: Copyright © DuMux Project contributors, see AUTHORS.md in root folder
-% SPDX-License-Identifier: CC-BY-4.0
-
-\section{Temporal Discretization and Solution Strategies}
-In this section, the temporal discretization as well as solution strategies (monolithic/sequential) are presented.
-
-\subsection{Temporal discretization}
-
-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
-\begin{equation}\label{eq:generalbalance}
-\frac{\partial m(u)}{\partial t} + \nabla\cdot\mathbf{f}(u, \nabla u) + q(u) = 0,
-\end{equation}
-seeking an unknown quantity $u$ in terms of storage $m$, flux $\mathbf{f}$ and source $q$.
-All available Dumux models can be written mathematically in form of \eqref{eq:generalbalance}
-with possibly vector-valued quantities $u$, $m$, $q$ and a tensor-valued flux $\mathbf{f}$.
-For the sake of simplicity, we assume scalar quantities $u$, $m$, $q$ and a vector-valued
-flux $\mathbf{f}$ in the notation below.
-
-For discretizing \eqref{eq:generalbalance}, 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
-\begin{equation}\label{eq:euler}
-\frac{\partial m(u_{k/k+1})}{\partial t}
-\approx \frac{m(u_{k+1}) - m(u_k)}{\Delta t_{k+1}}
-\end{equation}
-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.
-In case of the former, inserting \eqref{eq:euler} in \eqref{eq:generalbalance}
-at time $t_k$ leads to
-\begin{equation}\label{eq:expliciteuler}
-\frac{m(u_{k+1}) - m(u_k)}{\Delta t_{k+1}} + \nabla\cdot\mathbf{f}(u_k, \nabla u_k) + q(u_k) = 0,
-\end{equation}
-whereas the implicit Euler method is described as
-\begin{equation}\label{eq:impliciteuler}
-\frac{m(u_{k+1}) - m(u_k)}{\Delta t_{k+1}}
-+ \nabla\cdot\mathbf{f}(u_{k+1}, \nabla u_{k+1}) + q(u_{k+1}) = 0.
-\end{equation}
-Once the solution $u_k$ at time $t_k$ is known, it is straightforward
-to determine $m(u_{k+1})$ from \eqref{eq:expliciteuler},
-while attempting to do the same based on \eqref{eq:impliciteuler}
-involves the solution of a system of equations.
-On the other hand, the explicit method \eqref{eq:expliciteuler} is stable only
-if the time step size $\Delta t_{k+1}$ is below a certain limit that depends
-on the specific balance equation, whereas the implicit method \eqref{eq:impliciteuler}
-is unconditionally stable.
-
-\subsection{Solution strategies to solve equations}
-The governing equations of each model can be solved monolithically or sequentially.
-The basic idea of the sequential algorithm is to reformulate the
-equations of multi-phase flow into one equation for
-pressure and equations for phase/component/... transport. The pressure equation
-is the sum of the mass balance equations and thus considers the total flow of the
-fluid system. The new set of equations is considered as decoupled (or weakly coupled)
-and can thus be solved sequentially. The most popular sequential model is the
-fractional flow formulation for two-phase flow which is usually implemented applying
-an IMplicit Pressure Explicit Saturation algorithm (IMPES).
-In comparison to solving the equations monolithically, the sequential structure allows the use of
-different discretization methods for the different equations. The standard method
-used in the sequential algorithm is a cell-centered finite volume method. Further schemes,
-so far only available for the two-phase pressure equation, are cell-centered finite
-volumes with multi-point flux approximation (Mpfa-O method) and mimetic finite differences.
-An $h$-adaptive implementation of both sequential algorithms is provided for two dimensions.
--
GitLab