From: bangerth Date: Tue, 23 Nov 2010 02:18:43 +0000 (+0000) Subject: Add. X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=7ef874da8ff5057acc3c6c7353a75cf0939ab53e;p=dealii-svn.git Add. git-svn-id: https://svn.dealii.org/trunk@22843 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/examples/step-43/Makefile b/deal.II/examples/step-43/Makefile new file mode 100644 index 0000000000..a15735c1a2 --- /dev/null +++ b/deal.II/examples/step-43/Makefile @@ -0,0 +1,170 @@ +# $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========= $( $@" + @$(CXX) $(CXXFLAGS.g) -c $< -o $@ +./%.$(OBJEXT) : + @echo "==============optimized===== $( $@" + @$(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 diff --git a/deal.II/examples/step-43/doc/builds-on b/deal.II/examples/step-43/doc/builds-on new file mode 100644 index 0000000000..b27e0d398f --- /dev/null +++ b/deal.II/examples/step-43/doc/builds-on @@ -0,0 +1 @@ +step-21 step-31 diff --git a/deal.II/examples/step-43/doc/intro.dox b/deal.II/examples/step-43/doc/intro.dox new file mode 100644 index 0000000000..70537a872f --- /dev/null +++ b/deal.II/examples/step-43/doc/intro.dox @@ -0,0 +1,547 @@ +
+ + +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. + + + +

Introduction

+ +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: +
    +
  • Higher order spatial discretizations +
  • Adaptive mesh refinement +
  • Adaptive time stepping method +
  • Operator splitting method +
  • Efficient solver and preconditioning method +
+ +

Advection-dominated two-phase flow mathematical model.

+ +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}$. + +

Adaptive operator splitting and time stepping.

+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_pTime discretization. +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} + +

Weak form, space discretization for the pressure-velocity part.

+ +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. + +

Stabilization, weak form and space discretization for the saturation transport equation.

+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_{(-)}$. + +

Adaptive mesh refinement.

+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. + +

Linear system and its preconditioning.

+ +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. + +

The test cases.

+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 Trilinos 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: + + + + + + + + + + + + +
PARAMETER SYMBOL VALUE UNIT
Porosity $\epsilon$ 1.0 -
Viscosity (wetting) $\mu_w$ 0.2 $kg \cdot m^{-1} \cdot sec^{-1}$
Viscosity (nonwetting) $\mu_{nw}$ 1.0 $kg \cdot m^{-1} \cdot sec^{-1}$
Stabilization exponent $\alpha$ 1.0 -
Stabilization constant $\beta$ 2D: 0.3; 3D: 0.27 -
Normalization constant $c_R$ 1.0 -
Number of high-permeability regions $N$ 50; 200 -
Operator splitting threshold $\theta^\ast$ 5.0 -
+ + +

List of references

+ + +
    +
  1. +CC Chueh, W Bangerth, and N Djilali. +
    An h-adaptive operator splitting method for two-phase flow in 3D + heterogeneous porous media. +
    Submitted to Communications in Computational Physics. + +
  2. +F Brezzi and M Fortin. +
    Mixed and Hybrid Finite Element Methods. +
    Springer-Verlag, 1991. + +
  3. +Z Chen. +
    Finite Element Methods and Their Applications. +
    Springer, 2005. + +
  4. +JL Guermond and R Pasquetti. +
    Entropy-based nonlinear viscosity for fourier approximations of + conservation laws. +
    Comptes Rendus Mathematique, 346(13-14):801-806, 2008. + +
  5. +CC Chueh, M Secanell, W Bangerth, and N Djilali. +
    Multi-level adaptive simulation of transient two-phase flow in + heterogeneous porous media. +
    Computers and Fluids, 39:1585-1596, 2010. + +
  6. +Y Saad and MH Schultz. +
    Gmres: A generalized minimal residual algorithm for solving + nonsymmetric linear systems. +
    SIAM Journal on Scientific and Statistical Computing, + 7(3):856-869, 1986. + +
  7. +F Zhang. +
    The Schur Complement and its Applications. +
    Springer, 2005. + +
  8. +D Silvester and A Wathen. +
    Fast iterative solution of stabilised Stokes systems part ii: Using + general block preconditioners. +
    SIAM Journal on Numerical Analysis, 31(5):1352-1367, 1994. + +
  9. +GH Golub and CF van Loan. +
    Matrix Computations. +
    3rd Edition, Johns Hopkins, 1996. + +
  10. +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. +
    An overview of the trilinos project. +
    ACM Trans. Math. Softw., 31(3):397-423, 2005. + +
  11. +SE Buckley and MC Leverett. +
    Mechanism of fluid displacements in sands. +
    AIME Trans., 146:107-116, 1942. + +
diff --git a/deal.II/examples/step-43/doc/kind b/deal.II/examples/step-43/doc/kind new file mode 100644 index 0000000000..e62f4e7222 --- /dev/null +++ b/deal.II/examples/step-43/doc/kind @@ -0,0 +1 @@ +fluids diff --git a/deal.II/examples/step-43/doc/results.dox b/deal.II/examples/step-43/doc/results.dox new file mode 100644 index 0000000000..b5eaba9377 --- /dev/null +++ b/deal.II/examples/step-43/doc/results.dox @@ -0,0 +1,2 @@ +

Results

+ diff --git a/deal.II/examples/step-43/doc/tooltip b/deal.II/examples/step-43/doc/tooltip new file mode 100644 index 0000000000..48a5223056 --- /dev/null +++ b/deal.II/examples/step-43/doc/tooltip @@ -0,0 +1 @@ +Efficient ways to solve two-phase flow problems on adaptive meshes in 2d and 3d. diff --git a/deal.II/examples/step-43/step-40.cc b/deal.II/examples/step-43/step-40.cc new file mode 100644 index 0000000000..3a2cfde492 --- /dev/null +++ b/deal.II/examples/step-43/step-40.cc @@ -0,0 +1,15 @@ +/* $Id$ */ +/* Author: Chih-Che Chueh, University of Victoria, 2010 */ +/* Wolfgang Bangerth, Texas A&M University, 2010 */ + +/* $Id$ */ +/* */ +/* Copyright (C) 2010 by Chih-Che Chueh and the deal.II authors */ +/* */ +/* This file is subject to QPL and may not be distributed */ +/* without copyright and license information. Please refer */ +/* to the file deal.II/doc/license.html for the text and */ +/* further information on this license. */ + +int main () +{}