--- /dev/null
+# $Id$
+
+
+# For the small projects Makefile, you basically need to fill in only
+# four fields.
+#
+# The first is the name of the application. It is assumed that the
+# application name is the same as the base file name of the single C++
+# file from which the application is generated.
+target = $(basename $(shell echo step-*.cc))
+
+# The second field determines whether you want to run your program in
+# debug or optimized mode. The latter is significantly faster, but no
+# run-time checking of parameters and internal states is performed, so
+# you should set this value to `on' while you develop your program,
+# and to `off' when running production computations.
+debug-mode = on
+
+
+# As third field, we need to give the path to the top-level deal.II
+# directory. You need to adjust this to your needs. Since this path is
+# probably the most often needed one in the Makefile internals, it is
+# designated by a single-character variable, since that can be
+# reference using $D only, i.e. without the parentheses that are
+# required for most other parameters, as e.g. in $(target).
+D = ../../
+
+
+# The fourth field specifies the names of data and other files that
+# shall be deleted when calling `make clean'. Object and backup files,
+# executables and the like are removed anyway. Here, we give a list of
+# files in the various output formats that deal.II supports.
+clean-up-files = *gmv *gnuplot *gpl *eps *pov *vtk *ucd *.d2
+
+
+
+# The last field specifies the name of the input file that passes the
+# parameters to the program.
+run-parameters = input.prm
+
+
+
+
+
+#
+#
+# Usually, you will not need to change anything beyond this point.
+#
+#
+# The next statement tells the `make' program where to find the
+# deal.II top level directory and to include the file with the global
+# settings
+include $D/common/Make.global_options
+
+################################################################
+# This example program will only work if Trilinos is installed. If this
+# is not the case, then simply redefine the main targets to do nothing
+ifneq ($(USE_CONTRIB_TRILINOS),yes)
+default run clean:
+ @echo
+ @echo "==========================================================="
+ @echo "= This program cannot be compiled without Trilinos. Make="
+ @echo "= sure you have Trilinos installed and detected during ="
+ @echo "= configuration of deal.II ="
+ @echo "==========================================================="
+ @echo
+
+else
+#
+################################################################
+
+
+
+
+# Since the whole project consists of only one file, we need not
+# consider difficult dependencies. We only have to declare the
+# libraries which we want to link to the object file. deal.II has two
+# libraries: one for the debug mode version of the
+# application and one for optimized mode.
+libs.g := $(lib-deal2.g)
+libs.o := $(lib-deal2.o)
+
+
+# We now use the variable defined above to switch between debug and
+# optimized mode to select the set of libraries to link with. Included
+# in the list of libraries is the name of the object file which we
+# will produce from the single C++ file. Note that by default we use
+# the extension .g.o for object files compiled in debug mode and .o for
+# object files in optimized mode (or whatever local default on your
+# system is instead of .o)
+ifeq ($(debug-mode),on)
+ libraries = $(target).g.$(OBJEXT) $(libs.g)
+else
+ libraries = $(target).$(OBJEXT) $(libs.o)
+endif
+
+
+# Now comes the first production rule: how to link the single object
+# file produced from the single C++ file into the executable. Since
+# this is the first rule in the Makefile, it is the one `make' selects
+# if you call it without arguments.
+$(target)$(EXEEXT) : $(libraries)
+ @echo ============================ Linking $@
+ @$(CXX) -o $@ $^ $(LIBS) $(LDFLAGS)
+
+
+# To make running the application somewhat independent of the actual
+# program name, we usually declare a rule `run' which simply runs the
+# program. You can then run it by typing `make run'. This is also
+# useful if you want to call the executable with arguments which do
+# not change frequently. You may then want to add them to the
+# following rule:
+run: $(target)$(EXEEXT)
+ @echo ============================ Running $<
+ @./$(target)$(EXEEXT) $(run-parameters)
+
+
+# As a last rule to the `make' program, we define what to do when
+# cleaning up a directory. This usually involves deleting object files
+# and other automatically created files such as the executable itself,
+# backup files, and data files. Since the latter are not usually quite
+# diverse, you needed to declare them at the top of this file.
+clean:
+ -rm -f *.$(OBJEXT) *~ Makefile.dep $(target)$(EXEEXT) $(clean-up-files)
+
+
+# Since we have not yet stated how to make an object file from a C++
+# file, we should do so now. Since the many flags passed to the
+# compiler are usually not of much interest, we suppress the actual
+# command line using the `at' sign in the first column of the rules
+# and write the string indicating what we do instead.
+./%.g.$(OBJEXT) :
+ @echo "==============debug========= $(<F) -> $@"
+ @$(CXX) $(CXXFLAGS.g) -c $< -o $@
+./%.$(OBJEXT) :
+ @echo "==============optimized===== $(<F) -> $@"
+ @$(CXX) $(CXXFLAGS.o) -c $< -o $@
+
+
+# The following statement tells make that the rules `run' and `clean'
+# are not expected to produce files of the same name as Makefile rules
+# usually do.
+.PHONY: run clean
+
+
+# Finally there is a rule which you normally need not care much about:
+# since the executable depends on some include files from the library,
+# besides the C++ application file of course, it is necessary to
+# re-generate the executable when one of the files it depends on has
+# changed. The following rule creates a dependency file
+# `Makefile.dep', which `make' uses to determine when to regenerate
+# the executable. This file is automagically remade whenever needed,
+# i.e. whenever one of the cc-/h-files changed. Make detects whether
+# to remake this file upon inclusion at the bottom of this file.
+#
+# If the creation of Makefile.dep fails, blow it away and fail
+Makefile.dep: $(target).cc Makefile \
+ $(shell echo $D/include/deal.II/*/*.h)
+ @echo ============================ Remaking $@
+ @$D/common/scripts/make_dependencies $(INCLUDE) -B. $(target).cc \
+ > $@ \
+ || (rm -f $@ ; false)
+ @if test -s $@ ; then : else rm $@ ; fi
+
+# To make the dependencies known to `make', we finally have to include
+# them:
+include Makefile.dep
+
+
+endif # CONTRIB_USE_TRILINOS
--- /dev/null
+<br>
+
+<i>
+This program was contributed by Chih-Che Chueh (University of Victoria) and
+Wolfgang Bangerth. Results from this program are used and discussed in the
+following publications:
+- Chih-Che Chueh, Marc Secanell, Wolfgang Bangerth, Ned Djilali. Multi-level
+ adaptive simulation of transient two-phase flow in heterogeneous porous
+ media. Computers & Fluids, 39:1585-1596, 2010
+- Chih-Che Chueh, Wolfgang Bangerth, Ned Djilali. An h-adaptive operator
+ splitting method for two-phase flow in 3D heterogeneous porous
+ media. Submitted to Communications in Computational Physics.
+
+The implementation discussed here uses and extends
+parts of the step-21, step-31 and step-33 tutorial programs.
+
+The work of the first author was funded through the Canada Research Chairs
+Program and the MITACS Network of Centres of Excellence. Parts of the work by
+the second author were funded through Award No. KUS-C1-016-04, made by the King
+Abdullah University of Science and Technology, and through an Alfred P. Sloan
+Research Fellowship.
+</i>
+
+
+<a name="Intro"></a> <h1>Introduction</h1>
+
+The simulation of multiphase flow in porous media is a ubiquitous problem, but it faces two major difficulties: numerical accuracy and efficiency. In this tutorial, in order to overcome these two problems, there are five areas which we are trying to improve for a high performance simulator:
+<ul>
+<li> Higher order spatial discretizations
+<li> Adaptive mesh refinement
+<li> Adaptive time stepping method
+<li> Operator splitting method
+<li> Efficient solver and preconditioning method
+</ul>
+
+<h3>Advection-dominated two-phase flow mathematical model.</h3>
+
+We consider the flow of a two-phase immiscible, incompressible
+fluid. Capillary and gravity effects are neglected, and that viscous
+effects are assumed dominant. The governing equations for such a
+flow that are identical to those used in step-21 are then
+@f{align*}
+ \mathbf{u}_t &= - \mathbf{K} \lambda_t \left(S\right) \nabla p, \\
+ \nabla \cdot \mathbf{u}_t &= q, \\
+ \epsilon \frac{\partial S}{\partial t} + \nabla \cdot \left( \mathbf{u}_t F\left( S \right) \right)&=0,
+@f}
+where $S$ is the saturation (volume fraction) of the second (wetting) phase, $p$ is the pressure, $\mathbf{K}$ is the permeability tensor, $\lambda_t$ is the total mobility, $\epsilon$ is the porosity, $F$ is the fractional flow of the wetting phase, $q$ is the source term and $\mathbf{u}_t$ is the total velocity. The total mobility, fractional flow of the wetting phase and total velocity are respectively given by
+@f{align*}
+ \lambda_t(S)&= \lambda_w + \lambda_{nw} = \frac{k_{rw}(S)}{\mu_w} + \frac{k_{rnw}(S)}{\mu_{nw}}, \\
+ F(S) &= \frac{\lambda_w}{\lambda_t} = \frac{\lambda_w}{\lambda_w + \lambda_{nw}} = \frac{k_{rw}(S)/\mu_w}{k_{rw}(S)/\mu_w + k_{rnw}(S)/\mu_{nw}}, \\
+ \mathbf{u}_t &= \mathbf{u}_w + \mathbf{u}_{nw} = -\lambda_t(S)\mathbf{K} \cdot \nabla p,
+@f}
+where subscripts $w, nw$ represent the wetting and non-wetting phases,
+respectively.
+
+For convenience, the
+porosity $\epsilon$ in the saturation equation, which can be considered a
+scaling factor for the time variable, is set to
+one. Following a commonly used prescription for the dependence of the relative
+permeabilities $k_{rw}$ and $k_{rnw}$ on saturation, we use
+@f{align*}
+ k_{rw} &= S^2, \qquad\qquad
+ k_{rnw} &= \left( 1-S \right)^2.
+@f}
+
+The porous media equations above are
+augmented by initial conditions for the saturation and boundary conditions for
+the pressure. Since saturation and the gradient of the pressure uniquely
+determine the velocity, no boundary conditions are necessary for the velocity.
+Since the flow equations do not contain time derivatives, initial conditions for the velocity and pressure
+variables are not required. The flow field separates the boundary into inflow or outflow
+parts. Specifically,
+@f[
+ \mathbf{\Gamma}_{in}(t) = \left\{\vec{x} \in \partial \Omega:\vec{n} \cdot \vec{u}_t<0\right\},
+@f]
+and we arrive at a complete model by also imposing boundary values for the
+saturation variable on the inflow boundary $\mathbf{\Gamma}_{in}$.
+
+<h3>Adaptive operator splitting and time stepping.</h3>
+Based on facts that an implicit pressure-velocity part is computationally
+costly and the pressure and velocity depend only weakly on saturation, we
+don't have to solve it every saturation time step. Therefore, operator
+splitting is employed. However, with operator splitting that is not enough, we
+still needs an objective criterion to tell us when the pressure-velocity part
+needs to be solved at every saturation time step to control the amount of
+splitting error in a reasonable frame. Here, we use an a poteriori criterion
+(which detailed derivations and descriptions can be found in [Chueh, Bangerth
+and Djilali 2010]):
+@f{align*}
+ \theta(n,n_p)
+ =
+ \max_{\kappa\in{\mathbb T}}
+ \left(
+ \left\|
+ \frac 1{\lambda_t\left(S^{(n-1)}\right)}
+ - \frac 1{\lambda_t\left(S^{(n_p-1)}\right)} \right\|_{L^\infty(\kappa)}
+ \left\|\|\mathbf{K}^{-1}\|_1\right\|_{L^\infty(\kappa)}
+ \right).
+@f}
+where superscripts in parentheses denote the number of the saturation time step at which any quantity is defined and $n_p<n$ represents the last step where actually we computed the pressure and velocity.
+
+In short, the algorithm allows us to perform a number of
+saturation time steps of length $\Delta t_c^{(n)}$ until the criterion above
+tells us to re-compute velocity and pressure
+variables, leading to a macro time step of length
+@f[
+ \Delta t_p^{(n)} = \sum_{i=n_p+1}^{n} \Delta t_c^{(i)}.
+@f]
+We choose the length of (micro) steps subject to the Courant-Friedrichs-Lewy
+(CFL) restriction according to the criterion
+@f[
+ \Delta t_c = \frac{\textrm{min}_{K}h_{K}}{7 \|\vec{u}_t\|_{L^{\infty}\left(\Omega\right)}},
+@f]
+which we have confirmed to be stable for the choice of finite element and time
+stepping scheme for the saturation equation discussed below ($h_K$ denotes the
+diameter of cell $K$).
+The result is a scheme where neither micro nor macro time
+steps are of uniform length, and both are chosen adaptively.
+
+<h3>Time discretization.</h3>
+Using this time discretization, we obtain the following set of equations for each time step:
+@f{align*}
+ \mathbf{u}^{(n)}_t + \lambda_t\left(S^{(n-1)}\right) \mathbf{K} \nabla p^{(n)} =0, \\
+ \nabla \cdot \mathbf{u}^{(n)}_t = q, \\
+ \epsilon \left( \frac{S^{(n-1)}-S^{(n)}}{\Delta t^{(n)}_c} \right) + \mathbf{u}^{(n)}_t \cdot \nabla F\left(S^{(n-1)}\right) + F\left(S^{(n-1)}\right) \nabla \cdot \mathbf{u}^{(n)}_t =0,
+@f}
+where $\Delta
+t^{(n)}_c=t^{(n)}_c-t^{(n-1)}_c$ is the length of the $n$-th time step.
+
+Using the fact that $\nabla \cdot \mathbf{u}_t = q$, the time discrete
+saturation equation becomes
+@f{align*}
+ &\epsilon \left( \frac{S^{(n)}-S^{(n-1)}}{\Delta t^{(n)}_c} \right) + \mathbf{u}^{(n)}_t \cdot \nabla F\left(S^{(n-1)}\right) + F\left(S^{(n-1)}\right)q=0.
+@f}
+
+<h3>Weak form, space discretization for the pressure-velocity part.</h3>
+
+By multiplying the equations defining the total velocity $\mathbf u_t^{(n)}$ and
+the equation that expresses its divergence in terms of source terms, with test
+functions $\vec{v}$ and $w$
+respectively and then integrating terms by parts as necessary, the weak form
+of the problem reads: Find $\vec u, p$ so that for all test functions
+$\vec{v}, w$ there holds
+@f{gather*}
+ \left( \left( \mathbf{K} \lambda_t\left(S^{(n-1)}\right) \right)^{-1} \mathbf{u}^{(n)}_t, \vec{v}\right)_{\Omega} - \left(p^{(n)}, \nabla \cdot \vec{v}\right)_{\Omega} = -\left(p^{(n)}, \vec{n} \cdot \vec{v} \right)_{\partial \Omega}, \\
+ - \left( \nabla \cdot \mathbf{u}^{(n)}_t,w\right)_{\Omega} = - \big(q,w\big)_{\Omega}.
+@f}
+Here, $\vec{n}$ represents the unit outward normal vector to $\partial \Omega$ and the pressure $p^{(n+1)}$ can be prescribed weakly on the boundary $\partial \Omega$.
+
+We use continuous finite elements to discretize the velocity and pressure
+equations. Specifically, we use mixed finite elements to ensure high order approximation
+for both vector (e.g. a fluid velocity) and scalar variables (e.g. pressure)
+simultaneously. For saddle point problems, it is well established that
+the so-called Babuska-Brezzi or Ladyzhenskaya-Babuska-Brezzi (LBB) conditions
+[Brezzi 1991, Chen 2005] need to be satisfied to ensure stability of
+the pressure-velocity system. These stability conditions are satisfied in the
+present work by using elements for velocity that are one order higher than for
+the pressure, i.e. $u_h \in Q^d_{p+1}$ and $p_h \in Q_p$, where $p=1$, $d$ is
+the space dimension, and $Q_s$ denotes the space of tensor product Lagrange
+polynomials of degree $s$ in each variable.
+
+<h3>Stabilization, weak form and space discretization for the saturation transport equation.</h3>
+The chosen $Q_1$ elements for the saturation equation do not lead to a stable
+discretization without upwinding or other kinds of stabilization, and spurious
+oscillations will appear in the numerical solution. Adding an artificial
+diffusion term is is one approach to eliminating these oscillations
+[Chen 2005]. On the other hand, adding too much diffusion smears sharp
+fronts in the solution and suffers from grid-orientation difficulties
+[Chen 2005]. To avoid these effects, we use the artificial diffusion
+term proposed by [Guermond and Pasquetti 2008] and already
+validated in [Chueh et al. 2010].
+
+This method modifies the (discrete) weak form of the saturation equation
+to read
+@f{align*}
+ \left(\epsilon \frac{\partial S_h}{\partial t},\sigma_h\right)
+ -
+ \left(\mathbf{u}_t F\left( S_h \right),
+ \nabla \sigma\right)
+ +
+ \left(\mathbf n \cdot \mathbf{u}_t \hat F\left( S_h \right),
+ \sigma\right)_{\partial\Omega}
+ +
+ (\nu(S_h) \nabla S_h, \nabla \sigma_h)
+ &=0
+ \qquad
+ \forall \sigma_h,
+@f}
+where $\nu$ is the artificial diffusion parameter and $\hat F$ is an
+appropriately chosen numerical flux on the boundary of the domain (we choose
+the obvious full upwind flux for this).
+
+Following [Guermond and Pasquetti 2008], we use
+the parameter as a piecewise
+constant function set on each cell $K$ with the diameter $h_{K}$ as
+@f[
+ \nu(S)|_{K} = \beta \| \mathbf{u}_t \|_{L^{\infty}(K)} \textrm{min} \left\{ h_{K},h^{\alpha}_{K} \frac{\|\textrm{Res}(S)\|_{L^{\infty}(K)}}{c(\vec{u}_t,S)} \right\}
+@f]
+where $\alpha$ is a stabilization exponent and $\beta$ is a dimensionless
+user-defined stabilization constant. Following [Guermond and Pasquetti 2008]
+as well as the implementation in step-31, the velocity and saturation global
+normalization constant, $c(\vec{u}_t,S)$, and the residual $\textrm{Res}(S)$
+are respectively given by
+@f[
+ c(\vec{u}_t,S) = c_R \|\vec{u}_t\|_{L^{\infty}(\Omega)} \textrm{var}(S) | \textrm{diam} (\Omega) |^{\alpha - 2}
+@f]
+and
+@f[
+ \textrm{Res}(S) = \left( \epsilon \frac{\partial S}{\partial t} + \vec{u}_t \cdot \nabla F(S) + F(S)q \right) \cdot S^{\alpha - 1}
+@f]
+where $c_R$ is a second dimensionless user-defined constant,
+$\textrm{diam}(\Omega)$ is the diameter of the domain and $\textrm{var}(S) =
+\textrm{max}_{\Omega} S - \textrm{min}_{\Omega} S$ is the range of the present
+saturation values in the entire computational domain $\Omega$.
+
+This stabilization scheme has a number of advantages over simpler schemes such
+as finite volume (or discontinuous Galerkin) methods or streamline upwind
+Petrov Galerkin (SUPG) discretizations. In particular, the artificial
+diffusion term acts primarily in the vicinity of discontinuities
+since the residual is small in areas where the saturation is smooth. It
+therefore provides for a higher degree of accuracy. On the other hand, it is
+nonlinear since $\nu$ depends on the saturation $S$. We avoid this difficulty
+by treating all nonlinear terms explicitly, which leads to the following
+fully discrete problem at time step $n$:
+@f{align*}
+ &\left( \epsilon S^{(n)},\sigma\right)_{\Omega} - \Delta t^{(n)}_c \Big(F\left(S^{(n-1)}\right)\mathbf{u}^{*}_t,\nabla\sigma\Big)_{\Omega} + \Delta t^{(n)}_c \Big(F\left(S^{(n-1)}\right)\left(\vec{n}\cdot\mathbf{u}^{*}_t\right),\sigma\Big)_{\partial\Omega} \nonumber \\
+ & \quad = \left( \epsilon S^{(n-1)},\sigma\right)_{\Omega} - \Delta t^{(n)}_c \bigg(\nu\left(S^{(n-1)}\right)\nabla S^{(n-1)},\nabla\sigma\bigg)_{\Omega} \nonumber \\
+ & \qquad + \Delta t^{(n)}_c \bigg(\vec{n}\cdot\nu\left(S^{(n-1)}\right)\nabla S^{(n-1)},\sigma\bigg)_{\partial\Omega}
+@f}
+where $\mathbf{u}_t^{*}$ is the velocity linearly extrapolated from
+$\vec{u}^{(n_p)}_t$ and $\vec{u}^{(n_{pp})}_t$ to the current time $t^{(n)}$ if $\theta<\theta^*$ while $\mathbf{u}_t^{*}$ is $\vec{u}^{(n_p)}_t$ if $\theta>\theta^*$.
+Consequently, the equation is linear in $S_h^{(n)}$ and all that is required
+is to solve with a mass matrix on the saturation space.
+
+Since the Dirichlet boundary conditions for saturation are only imposed on the
+inflow boundaries, the third term on the left hand side of the equation above
+needs to be split further into two parts:
+@f{align*}
+ &\Delta t^{(n)}_c \Big(F\left(S^{(n-1)}\right)\left(\vec{n}\cdot\mathbf{u}^{(n)}_t\right),\sigma\Big)_{\partial\Omega} \nonumber \\
+ &\qquad= \Delta t^{(n)}_c \Big(F\left(S^{(n-1)}_{(+)}\right)\left(\vec{n}\cdot\mathbf{u}^{(n)}_{t(+)}\right),\sigma\Big)_{\partial\Omega_{(+)}} + \Delta t^{(n)}_c \Big(F\left(S^{(n-1)}_{(-)}\right)\left(\vec{n}\cdot\mathbf{u}^{(n)}_{t(-)}\right),\sigma\Big)_{\partial\Omega_{(-)}}
+@f}
+where $\partial\Omega_{(-)} = \left\{\vec{x} \in \partial\Omega : \vec{n}
+ \cdot \vec{u}_t<0\right\}$ and
+$\partial\Omega_{(+)} = \left\{\vec{x} \in \partial\Omega : \vec{n} \cdot
+ \vec{u}_t>0\right\}$ represent inflow and outflow boundaries,
+respectively. We choose values using an
+upwind formulation, i.e. $S^{(n-1)}_{(+)}$ and $\mathbf{u}^{(n)}_{t(+)}$
+correspond to the values taken from the present cell, while the values of
+$S^{(n-1)}_{(-)}$ and $\mathbf{u}^{(n)}_{t(-)}$ are those taken from the
+neighboring boundary $\partial\Omega_{(-)}$.
+
+<h3>Adaptive mesh refinement.</h3>
+choosing meshes adaptively to resolve sharp
+saturation fronts is an essential ingredient to achieve efficiency in our
+algorithm. Here, we use the same shock-type refinement approach used in
+[Chueh et al. 2010] to select those cells that should be refined or
+coarsened. The refinement indicator for each cell $K$ of the triangulation is
+computed by
+@f[
+ \eta_{K} = |\nabla S_h(\vec x_K)|
+@f]
+where $S_h(\vec x_K)$ is the discrete saturation variable evaluated at the
+center of cell $K$. This approach is analogous to ones frequently used in
+compressible flow problems, where density gradients are used to indicate
+refinement.
+
+<h3>Linear system and its preconditioning.</h3>
+
+Following the discretization of the governing equations
+discussed above, we
+obtain a linear system of equations in time step $(n)$ of the following form:
+@f[
+ \left(
+ \begin{array}{ccc}
+ \mathbf{M}^{\vec{u}} & \mathbf{B}^{T} & \mathbf{0} \\
+ \mathbf{B} & \mathbf{0} & \mathbf{0} \\
+ \mathbf{H} & \mathbf{0} & \mathbf{M}^{S}
+ \end{array}
+ \right)
+ \left(
+ \begin{array}{c}
+ \mathbf{U}^{(n)} \\
+ \mathbf{P}^{(n)} \\
+ \mathbf{S}^{(n)}
+ \end{array}
+ \right)
+ =
+ \left(
+ \begin{array}{c}
+ 0 \\
+ \mathbf{F}_{2} \\
+ \mathbf{F}_{3}
+ \end{array}
+ \right)
+@f]
+where the individual matrices and vectors are defined as follows using shape functions $\mathbf{v}_i$ for velocity, and $\phi_i$ for both pressure and saturation:
+@f{align*}
+ \mathbf{M}^{\vec{u}}_{ij}
+ &= \left( \left( \mathbf{K} \lambda_t\left(S^{(n-1)}\right) \right)^{-1}
+ \mathbf{v}_{i},\mathbf{v}_{j}\right)_{\Omega},
+ &
+ \mathbf{M}^{S}_{ij} &= \left(\epsilon \phi_i,\phi_j\right)_{\Omega}
+ \\
+ \mathbf{B}_{ij}
+ &= - \left( \nabla \cdot \mathbf{v}_{j},\phi_{i}\right)_{\Omega},
+ &
+ \mathbf{H}_{ij}
+ &= - \Delta t^{(n)}_c \Big( F\left(S^{(n-1)}\right) \mathbf{v}_i,\nabla\phi_j\Big)_{\Omega}
+ \\
+ \left(\mathbf{F}_{2}\right)_i
+ &= - \big(F\left(S^{(n-1)}\right)q,\phi_i\big)_{\Omega},
+@f}
+and $\mathbf{F}_{3}$ as given in the definition of the stabilized transport
+equation.
+
+The linear system above is of block triangular form if we consider the top
+left $2\times 2$ panel of matrices as one block. We can therefore first solve
+for the velocity and pressure (unless we decide to use $\mathbf U^{(n_p)}$ in
+place of the velocity)
+followed by a solve for the saturation variable. The first of these steps
+requires us to solve
+@f[
+ \left(
+ \begin{array}{cc}
+ \mathbf{M}^{\vec{u}} & \mathbf{B}^{T} \\
+ \mathbf{B} & \mathbf{0}
+ \end{array}
+ \right)
+ \left(
+ \begin{array}{c}
+ \mathbf{U}^{(n)} \\
+ \mathbf{P}^{(n)}
+ \end{array}
+ \right)
+ =
+ \left(
+ \begin{array}{c}
+ 0 \\
+ \mathbf{F}_{2}
+ \end{array}
+ \right)
+@f]
+We apply the Generalized Minimal Residual (GMRES) method [Saad and Schultz
+1986] to this linear system. The ideal preconditioner for the
+velocity-pressure system is
+@f{align*}
+\mathbf{P} =
+ \left(
+ \begin{array}{cc}
+ \mathbf{M}^{\vec{u}} & \mathbf{0} \\
+ \mathbf{B} & -\mathbf{S}
+ \end{array}
+ \right),
+ & \qquad
+ \mathbf{P}^{-1} =
+ \left(
+ \begin{array}{cc}
+ \left(\mathbf{M}^{\vec{u}}\right)^{-1} & \mathbf{0} \\
+ \mathbf{S}^{-1} \mathbf{B} \left(\mathbf{M}^{\vec{u}}\right)^{-1} & -\mathbf{S}^{-1}
+ \end{array}
+ \right)
+ @f}
+where
+$\mathbf{S}=\mathbf{B}\left(\mathbf{M}^{\vec{u}}\right)^{-1}\mathbf{B}^T$ is
+the Schur complement [Zhang 2005] of the system. This preconditioner is
+optimal since
+@f{align*}
+ \mathbf{P}^{-1}
+ \left(
+ \begin{array}{cc}
+ \mathbf{M}^{\vec{u}} & \mathbf{B}^{T} \\
+ \mathbf{B} & \mathbf{0}
+ \end{array}
+ \right)
+ =
+ \left(
+ \begin{array}{cc}
+ \mathbf{I} & \left(\mathbf{M}^{\vec{u}}\right)^{-1} \mathbf{B}^{T} \\
+ \mathbf{0} & \mathbf{I}
+ \end{array}
+ \right),
+@f}
+and consequently all eigenvalues are equal to one. GMRES with this
+preconditioner would then converge in one iteration.
+
+However, we cannot of course expect to use exact inverses of the
+velocity mass matrix and the Schur complement. We therefore follow the
+approach by [Silvester and Wathen 1994] originally proposed for
+the Stokes system. Adapting it to the current set of equations yield the
+preconditioner
+@f{align*}
+ \mathbf{\tilde{P}}^{-1} =
+ \left(
+ \begin{array}{cc}
+ \widetilde{\left(\mathbf{{M}}^{\vec{u}}\right)^{-1}}
+ & \mathbf{0} \\
+ \widetilde{\mathbf{{S}}^{-1}} \mathbf{B} \widetilde{\left(\mathbf{{M}}^{\vec{u}}\right)^{-1}} & -\widetilde{\mathbf{{S}}^{-1}}
+ \end{array}
+ \right)
+@f}
+where a tilde indicates an approximation of the exact inverse matrix. In
+particular, since $\left(\mathbf{{M}}^{\vec{u}}\right)^{-1}=\left( \left(
+ \mathbf{K} \lambda_t \right)^{-1}
+ \mathbf{v}_{i},\mathbf{v}_{j}\right)_{\Omega}$
+is a sparse symmetric and positive definite matrix, we choose for
+$\widetilde{\left(\mathbf{{M}}^{\vec{u}}\right)^{-1}}$ a single application of
+a sparse incomplete Cholesky decomposition of this matrix
+[Golub and Van Loan 1996].
+We note that the Schur complement that corresponds to the porous
+media flow operator in non-mixed form, $-\nabla \cdot [\mathbf K
+\lambda_t(S)]\nabla$ and
+$\mathbf{\tilde {S}} = \left( \left( \mathbf{K} \lambda_t \right) \nabla \phi_{i},\nabla \phi_{j}\right)_{\Omega}$
+should be a good approximation of the actual Schur complement matrix $\mathbf
+S$. Since both of these matrices are again symmetric and positive definite, we
+use an incomplete Cholesky decomposition of $\mathbf{\tilde S}$ for $\widetilde
+{\mathbf{{S}}^{-1}}$. It is important to note that $\mathbf{\tilde S}$ needs
+to be built with Dirichlet boundary conditions to ensure its invertibility.
+
+Once the velocity is available $\mathbf{U}^{(n)} \equiv \mathbf{u}^*_t$ (see
+[Chueh et al. 2010]), we can assemble $\mathbf{H}$ and
+$\mathbf{F}_{3}$ and solve for the saturations using
+@f{align*}
+ \mathbf{M}^{S} \mathbf{S}^{(n)} = \mathbf{F}_{3} - \mathbf{H} \mathbf{U}^{(n)}.
+@f}
+where the mass matrix $\mathbf{M}^{S}$ is solved by the conjugate gradient
+method, using an incomplete Cholesky decomposition as preconditioner once
+more.
+
+<h3>The test cases.</h3>
+The implementation discussed here uses and extends
+parts of the step-21, step-31 and step-33 tutorial programs of this library. We
+use the implementation of the incomplete Cholesky decomposition provided by
+the <code>Trilinos</code> library (see [Trilinos 2004]).
+
+We show numerical results that illustrate the
+efficiency and accuracy of our combined methods in solving the two-phase flow equations
+augmented by
+appropriate initial and boundary in conjunction with two different choices of the
+permeability model. In the problems considered, there is no internal source term ($q=0$).
+
+For simplicity, we choose $\Omega=[0,1]^d,d=2,3$, though all methods (as well
+as our implementation) should work equally well on general unstructured meshes.
+
+Initial conditions are only required for the saturation variable, and we
+choose $S(\mathbf{x},0)=0$, i.e. the porous medium is initially filled by the
+non-wetting phase. We prescribe a linear pressure on the boundaries:
+@f[
+ p(\vec{x},t) = 1 - x \qquad
+ \textrm{on} \quad \partial \Omega \times [0,T].
+@f]
+Pressure and saturation uniquely
+determine a velocity, and the velocity determines whether a boundary segment
+is an inflow or outflow boundary. On the inflow part of the boundary,
+$\mathbf{\Gamma}_{in}(t)$, we impose
+@f{align*}
+ S(\vec{x},t) = 1 \qquad & \textrm{on} \quad \mathbf{\Gamma}_{in}(t) \cap \left\{x = 0\right\}, \\
+ S(\vec{x},t) = 0 \qquad & \textrm{on} \quad \mathbf{\Gamma}_{in}(t) \backslash \left\{x = 0\right\}.
+@f}
+In other words, the domain is flooded by the wetting phase from the left.
+No boundary conditions for the saturation are required for the outflow parts
+of the boundary.
+
+All the numerical and physical parameters used for the 2D/3D
+cases are listed in the following table:
+
+<table align="center" class="tutorial" width="50%">
+<tr>
+ <td>PARAMETER </td><td>SYMBOL </td><td>VALUE </td><td>UNIT </td></tr><tr>
+ <td>Porosity </td><td>$\epsilon$ </td><td>1.0 </td><td>- </td></tr><tr>
+ <td>Viscosity (wetting) </td><td>$\mu_w$ </td><td>0.2 </td><td>$kg \cdot m^{-1} \cdot sec^{-1}$ </td></tr><tr>
+ <td>Viscosity (nonwetting) </td><td>$\mu_{nw}$ </td><td>1.0 </td><td>$kg \cdot m^{-1} \cdot sec^{-1}$ </td></tr><tr>
+ <td>Stabilization exponent </td><td>$\alpha$ </td><td>1.0 </td><td>- </td></tr><tr>
+ <td>Stabilization constant </td><td>$\beta$ </td><td>2D: 0.3; 3D: 0.27 </td><td>- </td></tr><tr>
+ <td>Normalization constant </td><td>$c_R$ </td><td>1.0 </td><td>- </td></tr><tr>
+ <td>Number of high-permeability regions </td><td>$N$ </td><td>50; 200 </td><td>- </td></tr><tr>
+ <td>Operator splitting threshold </td><td>$\theta^\ast$ </td><td>5.0 </td><td>- </td></tr>
+</table>
+
+
+<h3>List of references</h3>
+
+
+<ol>
+<li>
+CC Chueh, W Bangerth, and N Djilali.
+<br> An h-adaptive operator splitting method for two-phase flow in 3D
+ heterogeneous porous media.
+<br> <i>Submitted to Communications in Computational Physics</i>.
+
+<li>
+F Brezzi and M Fortin.
+<br> <i>Mixed and Hybrid Finite Element Methods</i>.
+<br> Springer-Verlag, 1991.
+
+<li>
+Z Chen.
+<br> <i>Finite Element Methods and Their Applications</i>.
+<br> Springer, 2005.
+
+<li>
+JL Guermond and R Pasquetti.
+<br> Entropy-based nonlinear viscosity for fourier approximations of
+ conservation laws.
+<br> <i>Comptes Rendus Mathematique</i>, 346(13-14):801-806, 2008.
+
+<li>
+CC Chueh, M Secanell, W Bangerth, and N Djilali.
+<br> Multi-level adaptive simulation of transient two-phase flow in
+ heterogeneous porous media.
+<br> <i>Computers and Fluids</i>, 39:1585-1596, 2010.
+
+<li>
+Y Saad and MH Schultz.
+<br> Gmres: A generalized minimal residual algorithm for solving
+ nonsymmetric linear systems.
+<br> <i>SIAM Journal on Scientific and Statistical Computing</i>,
+ 7(3):856-869, 1986.
+
+<li>
+F Zhang.
+<br> <i>The Schur Complement and its Applications</i>.
+<br> Springer, 2005.
+
+<li>
+D Silvester and A Wathen.
+<br> Fast iterative solution of stabilised Stokes systems part ii: Using
+ general block preconditioners.
+<br> <i>SIAM Journal on Numerical Analysis</i>, 31(5):1352-1367, 1994.
+
+<li>
+GH Golub and CF van Loan.
+<br> <i>Matrix Computations</i>.
+<br> 3rd Edition, Johns Hopkins, 1996.
+
+<li>
+MA Heroux, RA Bartlett, VE Howle, RJ Hoekstra, JJ Hu, TG Kolda, RB Lehoucq,
+ KR Long, RP Pawlowski, ET Phipps, AG Salinger, HK Thornquist, RS Tuminaro,
+ JM Willenbring, A Williams, and KS Stanley.
+<br> An overview of the trilinos project.
+<br> <i>ACM Trans. Math. Softw.</i>, 31(3):397-423, 2005.
+
+<li>
+SE Buckley and MC Leverett.
+<br> Mechanism of fluid displacements in sands.
+<br> <i>AIME Trans.</i>, 146:107-116, 1942.
+
+</ol>
projects / dealii-svn.git / commitdiff