Here, we create an executable called \texttt{my\_test} from a source file \texttt{mainfile.cc}.
The name of the test will also be \texttt{my\_test} (has to be unique). The last argument specifies a command - here, we just run the executbable \texttt{my\_test} with an input file \texttt{params.input}. For more advanced uses of
The name of the test will also be \texttt{my\_test} (has to be unique). The last argument specifies a command - here, we just run the executable \texttt{my\_test} with an input file \texttt{params.input}. For more advanced uses of
the \texttt{add\_dune\_test} macro, have a look at the \texttt{test} directory. A complete documentation is given under \url{https://www.dune-project.org/sphinx/core-2.5/}.
@@ -25,17 +25,17 @@ Gmsh is an open-source flexible grid generator for unstructured finite-element m
the provided GUI, we refer to the Gmsh documentation (\url{http://geuz.org/gmsh/doc/texinfo/gmsh.html}).
The MSH format can contain element and boundary markers defined on the grid. Thus, boundaries can be easily marked as, e.g., inflow boundaries
using Gmsh. Further, the format supports higher order elements. They can be used to create boundary parameterization supported by, e.g., the grid
using Gmsh. Further, the format supports higher order elements. They can be used to create boundary parametrization supported by, e.g., the grid
manager \texttt{UGGrid}.
An example can be found in \texttt{dumux/test\allowbreak/io/gridmanager}.
\subsubsection{Eclipse Grid Format}
The Eclipse Grid Format (GRDECL) is commonly used for corner-point grids.
The Eclipse Grid Format (GRDECL) is commonly used for corner-point grids.
Such grids consist of hexahedra, which are described by eight points on so-called pillars. A special feature of corner-point geometries is that points on pillars can degenerate, meaning that two neighboring points on a pillar can coincide. Furthermore, faces are, in general, bi-linear and cells can be non-convex. This allows for the accurate description of faults, layers, or wells, occurring in geological environments.
Furthermore, petrophysical properties can be defined (for each cell), by using eclipse-specific keywords, e.g. \texttt{PORO}, \texttt{PERMX}, \texttt{PERMY}.
Furthermore, petrophysical properties can be defined (for each cell), by using eclipse-specific keywords, e.g. \texttt{PORO}, \texttt{PERMX}, \texttt{PERMY}.
\Dumux supports the Eclipse Grid Format by using the \texttt{opm-grid} module (see (\url{https://opm-project.org}).
\Dumux supports the Eclipse Grid Format by using the \texttt{opm-grid} module (see (\url{https://opm-project.org}).
An example can be found in \texttt{dumux/test\allowbreak/porousmediumflow/2p/implicit/cornerpoint}.
Here the basic definitions, the general models concept, and a list of
models available in \Dumux are given. The actual differential equations
can be found in the localresiduals (see doxygen documentation of the
can be found in the localresiduals (see Doxygen documentation of the
model's \texttt{LocalResidual} class).
\subsection{Basic Definitions and Assumptions}
The basic definitions and assumptions are made, using the example
of a three-phase three-component system water-NAPL-gas
\cite{A3:class:2002a}. The modification for other multicomponent
\cite{A3:class:2002a}. The modification for other multi-component
systems is straightforward and can be found, e.\ g., in
\cite{A3:bielinski:2006,A3:acosta:2006}.
...
...
@@ -139,10 +139,10 @@ Here $p_i$ refers to the partial pressure of component i.
As an example, if two equal volumes of gas A and gas B are mixed, the volume of the mixture stays the same but the pressures add up (see Figure \ref{fig:dalton1}).
The density of the mixture, $\varrho$, can be calculated as follows:
...
...
@@ -165,10 +165,10 @@ V = \sum_{i}^{}V_i.
As an example, if two volumes of gas A and B at equal pressure are mixed, the pressure of the mixture stays the same, but the volumes add up (see Figure \ref{fig:dalton2}).
For standard Galerkin schemes, the weighting functions $W_j$ are chosen the same as the ansatz functions $N_j$. However, this does not yield a locally mass-conservative scheme.
Therefore, for the Box method, the weighting functions $W_j$ are chosen as
the piecewise constant functions over a
For standard Galerkin schemes, the weighting functions $W_j$ are chosen the same as the ansatz functions $N_j$. However, this does not yield a locally mass-conservative scheme.
Therefore, for the Box method, the weighting functions $W_j$ are chosen as
the piece-wise constant functions over a
control volume box $B_j$, i.e.
\begin{equation}
W_j(x) = \begin{cases}
1 &x \in B_j \\
0 &x \notin B_j.\\
\end{cases}
\label{eq:weightingFunctions}
W_j(x) = \begin{cases}
1 &x \in B_j \\
0 &x \notin B_j.\\
\end{cases}
\label{eq:weightingFunctions}
\end{equation}
Thus, the Box method is a Petrov-Galerkin scheme, where the weighting functions do not belong to the same function space than the ansatz functions.
Inserting definition \eqref{eq:weightingFunctions} into equation \eqref{eq:weightedResidual} and using the \textsc{Green-Gaussian} integral theorem results in
+ \int_{\partial B_j} F(\tilde u^{n+1}) \cdot\mathbf n
\;\mathrm{d}\Gamma_{B_j} - Q_j^{n+1}\: = 0.
\end{equation}
Equation \eqref{eq:discfin} has to be fulfilled for each box $B_j$.
\subsection{Cell Centered Finite Volume Methods -- A Short Introduction}\label{cc}
Cell-centered finite volume methods use the elements of the grid as control volumes.
For each control volume the discrete values are determined at the element/control
volume center (not required to be the barycenters).
volume center (not required to be the barycenters).
We consider a domain $\Omega\subset\mathbb{R}^d$, $d \in\{2, 3\}$ with boundary $\Gamma=\partial\Omega$. Within this section, we consider the following elliptic problem
\begin{equation}
...
...
@@ -186,8 +186,8 @@ We consider a domain $\Omega \subset \mathbb{R}^d$, $d \in \{ 2, 3 \}$ with boun
\end{aligned}
\end{equation}
Here, $\mathbf{\Lambda}=\mathbf{\Lambda}(\mathbf{x}, \mathbf{u})$ is a symmetric and positive definite tensor of second rank (e.g. permeability, diffusivity, etc.), $u = u (\mathbf{x})$ is unknown and $q = q(\mathbf{x}, \mathbf{u})$ is a source/sink.
We denote by $\mathcal{M}$ the mesh that results from the division of the domain $\Omega$ into $n_e$ control volumes $K \subset\Omega$. Each $K$ is a polygonal open set such that $K \cap L =\emptyset, \forall{K \neq L}$ and $\overline{\Omega}=\cup_{K \in\mathcal{M}}\overline{K}$.
Here, $\mathbf{\Lambda}=\mathbf{\Lambda}(\mathbf{x}, \mathbf{u})$ is a symmetric and positive definite tensor of second rank (e.g. permeability, diffusivity, etc.), $u = u (\mathbf{x})$ is unknown and $q = q(\mathbf{x}, \mathbf{u})$ is a source/sink.
We denote by $\mathcal{M}$ the mesh that results from the division of the domain $\Omega$ into $n_e$ control volumes $K \subset\Omega$. Each $K$ is a polygonal open set such that $K \cap L =\emptyset, \forall{K \neq L}$ and $\overline{\Omega}=\cup_{K \in\mathcal{M}}\overline{K}$.
For the derivation of the finite-volume formulation we integrate the first equation of \eqref{eq:elliptic} over a control volume $K$ and apply the Gauss divergence theorem:
...
...
@@ -201,17 +201,17 @@ Splitting the control volume boundary $\partial K$ into a finite number of faces
where $F_{K, \sigma}$ is the discrete flux through face $\sigma$ flowing out of cell $K$ and $Q_K :=\int_K q \,\mathrm{d}x$ is the integrated source/sink term. Equation \eqref{eq:ccdisc} is the typical cell-centered finite-volume formulation.
Finite-volume schemes differ in the way how the term
$(\mathbf{\Lambda}_K \nabla u )\cdot\mathbf{n}$ is approximated (i.e. the choice of the fluxes $F_{K, \sigma}$). Using the symmetry of the tensor $\mathbf{\Lambda}_K$, this term can be rewritten as
$\nabla u \cdot\mathbf{\Lambda}_K\mathbf{n}$, which corresponds to the directional derivative of $u$ in co-normal direction $\mathbf{\Lambda}_K\mathbf{n}$.
where $F_{K, \sigma}$ is the discrete flux through face $\sigma$ flowing out of cell $K$ and $Q_K :=\int_K q \,\mathrm{d}x$ is the integrated source/sink term. Equation \eqref{eq:ccdisc} is the typical cell-centered finite-volume formulation.
Finite-volume schemes differ in the way how the term
$(\mathbf{\Lambda}_K \nabla u )\cdot\mathbf{n}$ is approximated (i.e. the choice of the fluxes $F_{K, \sigma}$). Using the symmetry of the tensor $\mathbf{\Lambda}_K$, this term can be rewritten as
$\nabla u \cdot\mathbf{\Lambda}_K\mathbf{n}$, which corresponds to the directional derivative of $u$ in co-normal direction $\mathbf{\Lambda}_K\mathbf{n}$.
In the following, the main ideas of the two-point flux approximation and the multi-point flux approximation methods are briefly described. Hereby, we restrict the discussion to the two-dimensional case.
Please also note that other types of equations, e.g. instationary parabolic problems, can be discretized by applying some time discretization scheme to the time derivatives and by using the finite-volume scheme for the flux discretization. For simplicity the discussion is restricted to the elliptic problem \eqref{eq:elliptic}.
\subsubsection{Tpfa Method}\label{cc_tpfa}
The linear two-point flux approximation is a simple but robust cell-centered finite-volume scheme, which is commonly used in commercial software.
This scheme can be derived by using the conormal decomposition, which reads
The linear two-point flux approximation is a simple but robust cell-centered finite-volume scheme, which is commonly used in commercial software.
This scheme can be derived by using the co-normal decomposition, which reads
@@ -226,7 +226,7 @@ with the tensor $\mathbf{\Lambda}_K$ associated with control volume $K$, the dis
\end{figure}
With these notations, it follows that for each cell $K$ and face $\sigma$
With these notations, it follows that for each cell $K$ and face $\sigma$
\begin{equation}
\nabla u \cdot\mathbf{\Lambda}_K \mathbf{n}_{K, \sigma} = t_{K,\sigma}\nabla u \cdot\mathbf{d}_{K,\sigma} + \nabla u \cdot\mathbf{d}^{\bot}_{K,\sigma}.
\end{equation}
...
...
@@ -235,7 +235,7 @@ For the Tpfa scheme, the second part in the above equation is neglected. By usin
By neglecting the orthogonal term, the consistency of the scheme is lost for general grids, where $\nabla u \cdot\mathbf{d}^{\bot}_{K,\sigma}\not=0$. The consistency is achieved only for so-called K-orthogonal grids for which $\mathbf{d}^{\bot}_{K,\sigma}=0$. For such grids we deduce that
By neglecting the orthogonal term, the consistency of the scheme is lost for general grids, where $\nabla u \cdot\mathbf{d}^{\bot}_{K,\sigma}\not=0$. The consistency is achieved only for so-called K-orthogonal grids for which $\mathbf{d}^{\bot}_{K,\sigma}=0$. For such grids we deduce that
with $\tau_{K,\sigma} :=\mathbf{n}_{K, \sigma}\mathbf{\Lambda}_K\mathbf{n}_{K, \sigma}, \tau_{L,\sigma} :=\mathbf{n}_{L, \sigma}\mathbf{\Lambda}_L\mathbf{n}_{L, \sigma}$, $d_{K,\sigma}:=\mathbf{n}_{K, \sigma}\cdot\mathbf{d}_{K, \sigma}$, and $d_{L,\sigma}:=\mathbf{n}_{L, \sigma}\cdot\mathbf{d}_{L, \sigma}$. This reduces, for the case of scalar permeability, to a distance weighted harmonic averaging of permeabilities.
\subsubsection{Mpfa Method}\label{cc_mpfa}
Expressions for the face fluxes $F_{K, \sigma}$ are obtained by introducing intermediate face unknowns $u_\sigma$ in addition to the cell unknowns $u_K$ and enforcing the physically motivated continuity of fluxes and continuity of the solution across the faces. For a face $\sigma$ between the two polygons $K$ and $L$ these conditions read:
...
...
@@ -284,8 +284,8 @@ with $t_{K,\sigma},t_{L,\sigma}$ as defined in equation \eqref{eq:conormalDecTpf
In the following, a multi-point flux approximation method (Mpfa-O method), which was introduced in \citet{A3:aavatsmark:2002}, is presented. The main difference to the Tpfa scheme is the fact that a consistent discrete gradient is constructed, i.e. the term $\nabla u \cdot\mathbf{d}^{\bot}_{K,\sigma}$ is not neglected.
For this scheme, a dual grid is created by connecting the barycenters of the cells with the barycenters of the faces ($d=2$) or the barycenters of the faces and edges ($d=3$). This divides each cell into sub-control volumes $K^v$. Analogously, each face is sub-divided into sub-control volume faces $\sigma^v$, see Figure \ref{pc:interactionRegion_mpfa}. We allow for piecewise constant $\mathbf{\Lambda}$ (denoted as $\mathbf{\Lambda}_K$ for each cell $K$) and construct discrete gradients $\nabla_\mathcal{D}^{K^v} u$ (per sub-control volume $K^v$).
In the following, we restrict our discussion to the two-dimensional setup that is shown in Figure \ref{pc:interactionRegion_mpfa}.
For this scheme, a dual grid is created by connecting the barycenters of the cells with the barycenters of the faces ($d=2$) or the barycenters of the faces and edges ($d=3$). This divides each cell into sub-control volumes $K^v$. Analogously, each face is sub-divided into sub-control volume faces $\sigma^v$, see Figure \ref{pc:interactionRegion_mpfa}. We allow for piecewise constant $\mathbf{\Lambda}$ (denoted as $\mathbf{\Lambda}_K$ for each cell $K$) and construct discrete gradients $\nabla_\mathcal{D}^{K^v} u$ (per sub-control volume $K^v$).
In the following, we restrict our discussion to the two-dimensional setup that is shown in Figure \ref{pc:interactionRegion_mpfa}.
Here, the discrete gradients are constructed to be consistent such that the following conditions hold:
\begin{equation}
\nabla_\mathcal{D}^{K^v} u \cdot (\mathbf{x}_{\sigma^v_1}- \mathbf{x}_{K}) = u_{\sigma^v_1} - u_K, \quad\nabla_\mathcal{D}^{K^v} u \cdot (\mathbf{x}_{\sigma^v_3}- \mathbf{x}_{K}) = u_{\sigma^v_3} - u_K.
...
...
@@ -296,7 +296,7 @@ Thus, a discrete gradient (for sub-control volume $K^v$) that fulfills these co
\begin{bmatrix}
u_{\sigma^v_1} - u_K \\
u_{\sigma^v_3} - u_K
\end{bmatrix}, \qquad\text{ with }\;\mathbb{D}_{K^v} :=
\end{bmatrix}, \qquad\text{ with }\;\mathbb{D}_{K^v} :=
where $\mathbf{u}$ contains the three cell unknowns $u_K,u_L,u_M$ and $\mathbf{u}_{\sigma}$ the three face unknowns $u_{\sigma^v_1}, u_{\sigma^v_2}, u_{\sigma^v_3}$.
where $\mathbf{u}$ contains the three cell unknowns $u_K,u_L,u_M$ and $\mathbf{u}_{\sigma}$ the three face unknowns $u_{\sigma^v_1}, u_{\sigma^v_2}, u_{\sigma^v_3}$.
Inserting these face unknowns into the flux expression \eqref{eq:discreteFluxRef} yields
