]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Add Yan Li's step-21 for the IMPES problem
authorbangerth <bangerth@0785d39b-7218-0410-832d-ea1e28bc413d>
Tue, 17 Oct 2006 04:24:27 +0000 (04:24 +0000)
committerbangerth <bangerth@0785d39b-7218-0410-832d-ea1e28bc413d>
Tue, 17 Oct 2006 04:24:27 +0000 (04:24 +0000)
git-svn-id: https://svn.dealii.org/trunk@14010 0785d39b-7218-0410-832d-ea1e28bc413d

12 files changed:
deal.II/examples/step-21/Makefile [new file with mode: 0644]
deal.II/examples/step-21/doc/intro.dox [new file with mode: 0644]
deal.II/examples/step-21/doc/perm.png [new file with mode: 0644]
deal.II/examples/step-21/doc/pr4meshes.png [new file with mode: 0644]
deal.II/examples/step-21/doc/results.dox [new file with mode: 0644]
deal.II/examples/step-21/doc/solp4meshes.png [new file with mode: 0644]
deal.II/examples/step-21/doc/sols4meshes.png [new file with mode: 0644]
deal.II/examples/step-21/doc/solu14meshes.png [new file with mode: 0644]
deal.II/examples/step-21/doc/solu24meshes.png [new file with mode: 0644]
deal.II/examples/step-21/doc/yanli_project_introduction.tex [new file with mode: 0644]
deal.II/examples/step-21/doc/yanli_project_result.tex [new file with mode: 0644]
deal.II/examples/step-21/step-21.cc [new file with mode: 0644]

diff --git a/deal.II/examples/step-21/Makefile b/deal.II/examples/step-21/Makefile
new file mode 100644 (file)
index 0000000..1972f54
--- /dev/null
@@ -0,0 +1,156 @@
+# $Id: Makefile,v 1.4 2006/02/10 17:53:05 wolf Exp $
+
+
+# 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 proj.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 last 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
+
+
+
+
+#
+#
+# Usually, you will not need to change something beyond this point.
+#
+#
+# The next statement tell 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
+
+
+# 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, and there need
+# to be two sets of libraries: one for the debug mode version of the
+# application and one for the optimized mode. Here we have selected
+# the versions for 2d. Note that the order in which the libraries are
+# given here is important and that your applications won't link
+# properly if they are given in another order.
+#
+# You may need to augment the lists of libraries when compiling your
+# program for other dimensions, or when using third party libraries
+libs.g   = $(lib-deal2-2d.g) \
+          $(lib-deal2-3d.g) \
+          $(lib-lac.g)      \
+           $(lib-base.g)
+libs.o   = $(lib-deal2-2d.o) \
+          $(lib-deal2-3d.o) \
+          $(lib-lac.o)      \
+           $(lib-base.o)
+
+
+# We now use the variable defined above which 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 the 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) : $(libraries)
+       @echo ============================ Linking $@
+       @$(CXX) -o $@$(EXEEXT) $^ $(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)
+       @echo ============================ Running $<
+       @./$(target)$(EXEEXT)
+
+
+# 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 to created 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/*/*.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
+
+
diff --git a/deal.II/examples/step-21/doc/intro.dox b/deal.II/examples/step-21/doc/intro.dox
new file mode 100644 (file)
index 0000000..f693a95
--- /dev/null
@@ -0,0 +1 @@
+<a name="Intro"></a> <h1>Introduction</h1>
diff --git a/deal.II/examples/step-21/doc/perm.png b/deal.II/examples/step-21/doc/perm.png
new file mode 100644 (file)
index 0000000..1c7fe71
Binary files /dev/null and b/deal.II/examples/step-21/doc/perm.png differ
diff --git a/deal.II/examples/step-21/doc/pr4meshes.png b/deal.II/examples/step-21/doc/pr4meshes.png
new file mode 100644 (file)
index 0000000..0f0afe5
Binary files /dev/null and b/deal.II/examples/step-21/doc/pr4meshes.png differ
diff --git a/deal.II/examples/step-21/doc/results.dox b/deal.II/examples/step-21/doc/results.dox
new file mode 100644 (file)
index 0000000..89faa40
--- /dev/null
@@ -0,0 +1 @@
+<a name="Results"></a> <h1>Results</h1>
diff --git a/deal.II/examples/step-21/doc/solp4meshes.png b/deal.II/examples/step-21/doc/solp4meshes.png
new file mode 100644 (file)
index 0000000..1a15adf
Binary files /dev/null and b/deal.II/examples/step-21/doc/solp4meshes.png differ
diff --git a/deal.II/examples/step-21/doc/sols4meshes.png b/deal.II/examples/step-21/doc/sols4meshes.png
new file mode 100644 (file)
index 0000000..9cce4b4
Binary files /dev/null and b/deal.II/examples/step-21/doc/sols4meshes.png differ
diff --git a/deal.II/examples/step-21/doc/solu14meshes.png b/deal.II/examples/step-21/doc/solu14meshes.png
new file mode 100644 (file)
index 0000000..8264339
Binary files /dev/null and b/deal.II/examples/step-21/doc/solu14meshes.png differ
diff --git a/deal.II/examples/step-21/doc/solu24meshes.png b/deal.II/examples/step-21/doc/solu24meshes.png
new file mode 100644 (file)
index 0000000..4b19717
Binary files /dev/null and b/deal.II/examples/step-21/doc/solu24meshes.png differ
diff --git a/deal.II/examples/step-21/doc/yanli_project_introduction.tex b/deal.II/examples/step-21/doc/yanli_project_introduction.tex
new file mode 100644 (file)
index 0000000..64779eb
--- /dev/null
@@ -0,0 +1,218 @@
+\documentclass[12pt]{article}
+\topmargin 0pt \oddsidemargin 0pt \evensidemargin 0pt
+\textwidth=14truecm \textheight=21.5truecm
+\renewcommand{\baselinestretch}{1.5}
+\begin{document}
+\author{Yan Li}
+\title{A Numerical Simulation for Two Phase Flow}
+\maketitle
+\section{Introduction}
+\subsection{Abstract}
+
+In this project, we propose a numerical simulation for two phase
+flow problem in porous media. The two phase flow system includes one
+elliptic equation and one nonlinear transport equation. We apply
+mixed finite element method and Discontinuous Galerkin method for
+this system. Some numerical results for two dimensional case are
+given by $RT_{0}\times DQ_{0}\times DQ_{0}$.
+\\
+The numerical computation is based on $dealII$. We use vector shape
+functions from step9, DG method from step12, mixed method and Schur
+complement from step20 and many many useful tools from the library.
+
+\subsection{Two Phase Flow Problem }
+The modeling of two phase flow in porous media is important for both
+environmental rededication and the management of petroleum
+reservoirs. Practical situations involving two phase flow include
+the dispersal of a nonaqueous phase liquid in an aquifer or the
+displacement of a non-aqueous heterogeneity on the flow and
+transport. Simulation models, if they are to provide realistic
+predictions, must accurately account for these effects.
+%However,
+%because permeability heterogeneity occurs at many different length
+%scales, numerical flow models cannot in general resolve all of the
+%scales of variation.Therefore, approaches are needed for
+%representing the effects of subgrid scale variations on larger scale
+%flow results. Typically, upscaled or multiscale models are employed
+%for such systems. \\
+In our project,we consider a kind of periodic permeability,our
+numerical result shows that the heterogeneity effects are simulated accurately.\\
+Consider two phase flow in a reservoir $\Omega$ under the assumption
+that the displacement is dominated by viscous effects; i.e. we
+neglect the effects of gravity, compressibility, and capillary
+pressure. Porosity will be considered to be constant. The two phase
+will be referred to as water and oil, designated by subscripts $w$
+and $o$, respectively. We write Darcy's for each phase as follows:
+\begin{eqnarray}
+\mathbf{u}_{j} = \frac{k_{rj}(S)}{\mu_{j}} \mathbf{K} \cdot \nabla p
+\end{eqnarray}
+\indent where, $\mathbf{u}_{j}$ is the phase velocity, $K$ is the
+permeability tensor, $k_{rj}$ is the relative permeability to phase
+$j$($j=o,w$),$S$ is the water saturation(volume fraction), $P$ is
+pressure and $\mu_{j}$ is the viscosity of phase $j$($j=o,w$).\\
+Combining Darcy's law with a statement of conservation of mass
+allows us to express the governing equations in terms of the
+so-called pressure and saturation equations:
+\begin{eqnarray}
+\nabla \cdot (\mathbf{K}(x,y)\lambda(S) \nabla p)= q(x,y) && \forall(x,y)\in\Omega\\
+ S_{t} + \mathbf{u} \cdot \nabla F(S) = 0&& \forall(x,y)\in\Omega
+\end {eqnarray}
+
+
+\indent where, $\lambda$ is the total mobility, f is the fractional
+flow of water, $q$ is a source term and $\mathbf{u}$ is the total
+velocity, which are respectively given by:
+$$\mathbf{u} =
+\mathbf{u}_{o} + \mathbf{u}_{w} = -\lambda(S) \mathbf{K}\cdot\nabla
+p$$
+$$\lambda(S) = \frac{k_{rw}(S)}{\mu_{w}}+\frac{k_{ro}(S)}{\mu_{o}}$$
+$$F(S) = \frac{k_{rw}(S)/\mu_{w}}{k_{rw}(S)/\mu_{w} + k_{ro}(S)/\mu_{o}}$$
+
+
+\subsection{Discretization}
+
+ For simplicity, in our project we will assume no
+source $q=0$ and the heterogeneous porous medium is isotropic
+$\mathbf{K}(x,y) =
+k(x,y) \mathbf{I}$. \\
+Our two dimensional numerical simulation will be done on unit cell
+$\Omega = [0,1]\times [0,1]$ for $t\in [0,T]$.
+\begin {eqnarray}
+\mathbf{u}(x,y)+\mathbf{K}(x,y)\lambda(S) \nabla p= 0 && \forall(x,y)\in\Omega, \forall t\in [0,T]\\
+\nabla \cdot\mathbf{u}(x,y)= 0 && \forall(x,y)\in\Omega, \forall t \in [0,T] \\
+S_{t} + \mathbf{u} \cdot \nabla F(S) = 0&& \forall(x,y)\in\Omega,
+\forall t \in [0,T]
+\end {eqnarray}
+ Boundary conditions are:
+\[
+\begin {array}{cr}
+p(x,y)=1 & \forall(x,y)\in \Gamma_{1}:=\{(x,y)\in \partial \Omega: x=0\}\\
+p(x,y)=0 & \forall(x,y)\in  \Gamma_{2}:=\{(x,y)\in \partial \Omega: x=1\}\\
+\mathbf{u}(x,y)\cdot \mathbf{n}=0 & \forall(x,y)\in
+\partial\Omega \setminus(\Gamma_{1}\bigcup \Gamma_{2})
+\end {array}
+\]
+\\
+Initial conditions are:
+\[
+\begin {array}{cr}
+S(x,y,t=0)= 1& \forall (x,y) \in \Gamma_{1}\\
+S(x,y,t=0) = 0 & \forall(x,y)\in \partial \Omega \setminus
+\Gamma_{1}
+\end {array}
+\]
+\\
+We apply mixed finite method on velocity and pressure. To be
+well-posed, we choose Raviart-Thomas spaces $RT_{k}$ for
+$\mathbf{u}$ and discontinuous elements of class $DQ_{k}$ for $p$,
+then the mixed
+system is:\\
+Find $(\mathbf{u},p)\in RT_{k}\times DQ_{k}$ such that:
+\begin {eqnarray}
+\sum_{\kappa}\{ \int _{\kappa}(K \lambda)^{-1} \mathbf{u}\cdot
+\mathbf{v} dx - \int_{\kappa} p \nabla \cdot \mathbf{v} dx\}
+ =- \int_{\Gamma _{1}} \mathbf{v}\cdot \mathbf{n}&&  \forall\mathbf{v}\in RT_{k}(\Omega)\\
+\sum_{\kappa}\{\int (\nabla \cdot \mathbf{u}) q dx\} = 0 && \forall
+q\in DQ_{k}(\Omega)
+\end {eqnarray}
+For saturation, we also use discontinuous finite element method.
+i.e. Find $S^{n+1} \in DQ_{k}$ such that for all $ \phi \in DQ_{k}$,
+the following formulation holds:
+\begin {eqnarray}
+\sum_{\kappa}\{\int_{\kappa}\frac{S^{n+1}-S^{n}}{\triangle t} \phi
+dx + \int_{\kappa} (\mathbf{u}^{n+1}\cdot \nabla F(S^{n})) \phi
+dx\} =0
+\end {eqnarray}
+Integrating by parts:
+\begin {eqnarray}
+\nonumber
+ \sum_{\kappa}\{\int_{\kappa}S^{n+1} \phi dx +\triangle t
+\int_{\partial \kappa}F(S^{n})( \mathbf{u}^{n+1}\cdot \mathbf{n} )
+\phi dx &-\triangle t\int_{\kappa}  F(S^{n})( \mathbf{u^{n+1}}\cdot
+\nabla
+\phi )dx\}\\
+&= \sum_{\kappa}\int_{\kappa} S^{n} \phi dx
+\end {eqnarray}
+
+\indent where,$\mathbf{n}$ denotes the unit outward normal to the
+boundary $\partial \kappa$. And here we can use $u^{n+1}$ instead of
+$u^{n}$ is because that we view $(u^{n+1},p^{n+1},S^{n+1})$ as
+a block vector,$u^{n+1}$ could be implement in the coefficient function for saturation.
+We believe the saturation is computed more accurately in this way.\\
+Considering the discontinuity of the discrete function $S_h$ on
+interelement faces, the flux $\mathbf{u}^{n+1}\cdot \mathbf{n} $ is
+computed as:
+ \begin{eqnarray}
+&&\int_{\partial \kappa}F(S^{n}) (\mathbf{u}^{n+1}\cdot \mathbf{n})
+\phi dx =\\
+\nonumber && \int_{\partial \kappa _{+}}
+F(S^{n,+})(\mathbf{u}^{n+1,+}\cdot \mathbf{n})\phi dx
++\int_{\partial \kappa _{-}} F(S^{n,-})(\mathbf{u}^{n+1,-}\cdot
+\mathbf{n})\phi dx
+\end{eqnarray}
+
+where, $\partial \kappa _{-}:= \{x\in
+\partial\kappa , \mathbf{u}(x) \cdot \mathbf{n}<0\}$ denotes the inflow boundary
+and$\partial \kappa _{+}:= \{\partial \kappa \setminus \partial
+\kappa_{-}\}$ is the outflow part of the boundary. By the
+discontinuity of$ S_{h}$ , $F(S^{n,-})$ takes the value of
+neighboring cell,$F(S^{n+})$ takes the value of cell $\kappa$.
+
+\subsection{Implementation}
+We use
+$dealII$ to implement our mixed and DG system. The main idea is same
+with step-20 but there are some new problems we have to consider:\\
+\indent $(1)$ We has the three blocks vector $(u,p,S)$ , in which
+all the functions are dependent on time. i.e. At each time step we
+need project the $solution$ into $old-solution$, using
+$old-solution$ to get a new $solution$.
+Keep doing this until the last time step;\\
+At time $t=t^{n+1}$ , suppose $old-solution=(u^{n},p^{n},S^{n})$ is
+known, in $assemble-system()$ part, we assemble system matrix as:
+\[
+\begin {array}{cccccccccccc}
+\lceil &M(S^{n}) &B^{T}& 0 &\rceil & \lceil& \mathbf{u}^{n+1}&\rceil& &\lceil& 0 &\rceil\\
+|      &B&    0 & 0 & |     &|      & p^{n+1} &|        &=&|     & q &|\\
+\lfloor&\triangle t \nabla F(S^n)&    0& I & \rfloor & \lfloor
+&S^{n+1} & \rfloor & & \lfloor& S^{n}& \rfloor
+\end {array}
+\]
+\\
+In $solve()$ part, we solve the first two equations independent of
+the third equation, since $M( S^n)$ is already known. As in step-20,
+using vector base functions, Schur complement with a
+preconditioner and CG method, we get $u^{n+1}$and $p^{n+1}$. \\
+Then, with the above $u^{n+1}$ and $p^{n+1}$, we could compute
+$S^{n+1}$ by :
+\begin {eqnarray}
+\sum_{\kappa}\int_{\kappa}S^{n+1} \phi dx&&=
+\sum_{\kappa}\{\int_{\kappa} S^{n} \phi dx+\Delta t\int_{\kappa}
+F(S^{n}) \mathbf{u^{n+1}}\cdot \nabla \phi dx\\
+\nonumber && -\Delta t \int_{\partial \kappa_{-}}F(S^{n,-})
+\mathbf{u}^{n+1,-}\cdot \mathbf{n} \phi dx -\Delta t \int_{\partial
+\kappa_{+}}F(S^{n}) \mathbf{u}^{n+1}\cdot \mathbf{n} \phi dx\}
+\end {eqnarray}
+Now, project solution $(u^{n+1},p^{n+1},S^{n+1})$ into
+$old-solution$, do the above process for next time step.\\
+ \indent
+$(2)$ The numerical flux term is related with neighbor cells.In our
+implementation $solve( )$, we do the following on each cell: \\
+For each face, compute the flux $\mathbf{u}\cdot F(S)$, the flux is
+negative means it is an in-flow face. Then if this in-flow face is
+on the boundary $\Gamma_{1}$:$F(S^{-})=F(1)$; If the in-flow
+face is not on boundary, $F(S^{-})=F(S|_{neighbor})$.
+Flux is positive means it is an out-flow face, we just use $ F(S)$ on current cell.\\
+All the other functions are commented in code, please see next part
+- the commented program.
+
+\subsection{Test Case}
+Our two phase flow are chosen as water and oil. The total mobility
+is : $$\lambda (S) = \frac{1.0}{\mu} S^2 +(1-S)^2$$ Permeability is
+:
+$$K(x,y)=\mathbf{k}(x,y)I=\frac{1.0}{2+1.99\sin(2\pi\frac{2x-y}{\epsilon})}
+I$$
+ Fractional flow of water is: $$F(S)=\frac{S^2}{S^2+\mu (1-S)^2}$$
+Choose $\epsilon=0.05$ , viscosity $\mu=0.2$. \\
+The resulting solution will be shown in result part.
+
+\end{document}
diff --git a/deal.II/examples/step-21/doc/yanli_project_result.tex b/deal.II/examples/step-21/doc/yanli_project_result.tex
new file mode 100644 (file)
index 0000000..699f6e6
--- /dev/null
@@ -0,0 +1,42 @@
+\documentclass[12pt]{article}
+\usepackage[dvips]{graphicx}
+
+\begin{document}
+\section{Result}
+
+For the test case in Introduction,we do numerical computation on four meshes:
+$$h_{1}=2^{-3}, h_{2}=2^{-4}, h_{3}=2^{-5}, h_{4}=2^{-6}$$
+
+Our permeability is $K(x,y)=\frac{1.0}{2+1.99*sin(2 \pi \frac{2x-y}{\epsilon})}$:\\
+
+\includegraphics[width=6 in]{perm.eps}
+
+Because of heterogeneity in the medium, the solution for velocity and saturatio are periodic.
+Compare the following pictures on four meshes at time $t=T$:\\
+Velocity-x:\\
+\includegraphics[width=5 in]{solu14meshes.eps}
+\\
+Velocity-y:\\
+\includegraphics[width=5 in]{solu24meshes.eps}
+\\
+Pressure:\\
+\includegraphics[width=5 in]{solp4meshes.eps}
+\\
+Saturation:\\
+\includegraphics[width=5 in]{sols4meshes.eps}
+\\
+From above comparation, we can see pressure is stable but velocity and saturation are not. 
+The reason is heterogeneity of the medium and some complexity fo the dynamic systems.
+By our direct numerical computation, only fine mesh solution is able to catch the subgrid properties.
+That means an accurate well-posed computation requires tremendous amount of computer memory and CPU time. 
+But usually ,it easily exceed the limit of today's computer resources.\\
+There are some alternative approaches have been developed. A common approach is to "scale up" a heterogeneous medium. 
+This method is to find an effective representation of permeablility on a coarse mesh so that the large scale flow can be correctly computed on this mesh.The computational cost is thus greatly reduced.
+\\
+At last,let's see Oil Production Rate on the boundary $\Gamma_{2}$:
+$$ PR(t)=1-\frac{\int_{\Gamma_{2}} (\mathbf{u}\cdot \mathbf{n})F(S)dx}{\int_{\Gamma_{2}} (\mathbf{u} \cdot \mathbf{n})dx}$$
+
+\includegraphics[width=6 in]{pr4meshes.eps}
+
+
+\end{document}
\ No newline at end of file
diff --git a/deal.II/examples/step-21/step-21.cc b/deal.II/examples/step-21/step-21.cc
new file mode 100644 (file)
index 0000000..ecb30a7
--- /dev/null
@@ -0,0 +1,1279 @@
+                                  // {Include files}
+
+                                // This program is an daptation of step-20
+                                // and includes some technique of DG method from step-12
+                                // We list include files in the order
+                                // base-lac-grid-dofs-fe-numerics.
+                                
+#include <base/quadrature_lib.h>
+#include <base/logstream.h>
+#include <base/function.h>
+#include <lac/block_vector.h>
+#include <lac/full_matrix.h>
+#include <lac/block_sparse_matrix.h>
+#include <lac/solver_cg.h>
+#include <lac/precondition.h>
+#include <grid/tria.h>
+#include <grid/grid_generator.h>
+#include <grid/tria_accessor.h>
+#include <grid/tria_iterator.h>
+#include <grid/grid_tools.h>
+#include <dofs/dof_handler.h>
+#include <dofs/dof_renumbering.h>
+#include <dofs/dof_accessor.h>
+#include <dofs/dof_tools.h>
+#include <dofs/dof_constraints.h>
+#include <fe/fe_q.h>
+                                  //The Discontinuous Galerkin finite element is declared: 
+#include <fe/fe_dgq.h>
+
+#include <fe/fe_system.h>
+#include <fe/fe_values.h>
+#include <fe/mapping_q1.h>
+#include <numerics/vectors.h>
+#include <numerics/matrices.h>
+#include <numerics/data_out.h>
+#include <fstream>
+#include <iostream>
+#include <sstream>
+                                 // The Raviart-Thomas finite element is declared:
+#include <fe/fe_raviart_thomas.h>
+
+                                // In this program, we use a tensorial
+                                // coefficient. Since it may have a
+                                // spatial dependence, we consider it
+                                // a tensor-valued function. The
+                                // following include file provides
+                                // the ``TensorFunction'' class that
+                                // offers such functionality:
+#include <base/tensor_function.h>
+
+
+                                 // {The ``TwoPhaseFlowProblem'' class template}
+                                 
+
+template <int dim>
+class TwoPhaseFlowProblem 
+{
+  public:
+    TwoPhaseFlowProblem (const unsigned int degree);
+    void run ();
+    
+  private:
+    void make_grid_and_dofs ();
+    void assemble_system ();
+    void solve ();
+    void compute_errors () const;
+    void output_results (const unsigned int timestep_number) const;
+
+    Vector<double> evaluate_solution (const Point<dim> &point) const;
+    
+    const unsigned int   degree;
+    
+    Triangulation<dim>   triangulation;
+    FESystem<dim>        fe;
+    DoFHandler<dim>      dof_handler;
+
+
+    BlockSparsityPattern      sparsity_pattern;
+    BlockSparseMatrix<double> system_matrix;
+
+    const unsigned int n_refinement_steps;
+    
+    double time_step;
+    double epsilon;
+    double vis;    
+    double vfs_out;
+    double v_out;
+    BlockVector<double>       solution;
+    BlockVector<double>       old_solution;
+    BlockVector<double>       system_rhs;
+    
+    
+};
+
+
+                                //{Right hand side, boundary values and initial values}
+                                
+                                  // we define the template for pressure right-hand side(source function)
+                                 //and boundary values for pressure and saturation
+                                 // initial values for saturation.
+
+template <int dim>
+class RightHandSide : public Function<dim> 
+{
+  public:
+    RightHandSide () : Function<dim>(1) {};
+    
+    virtual double value (const Point<dim>   &p,
+                         const unsigned int  component = 0) const;
+};
+
+
+
+template <int dim>
+class PressureBoundaryValues : public Function<dim> 
+{
+  public:
+    PressureBoundaryValues () : Function<dim>(1) {};
+    
+    virtual double value (const Point<dim>   &p,
+                         const unsigned int  component = 0) const;
+};
+
+
+template <int dim>
+class SaturationBoundaryValues : public Function<dim> 
+{
+  public:
+    SaturationBoundaryValues () : Function<dim>(dim+2) {};
+    
+    virtual void vector_value (const Point<dim> &p, 
+                              Vector<double>   &value) const;
+};
+
+
+template <int dim>
+class InitialValues : public Function<dim> 
+{
+  public:
+    InitialValues () : Function<dim>(dim+2) {};
+    
+    virtual double value (const Point<dim>   &p,
+                         const unsigned int  component = 0) const;
+
+    virtual void vector_value (const Point<dim> &p, 
+                              Vector<double>   &value) const;
+
+};
+
+
+
+
+                                // And then we also have to define
+                                // these respective functions, of
+                                // course. Given our discussion in
+                                // the introduction of how the
+                                // solution should look like, the
+                                // following computations should be
+                                // straightforward:
+template <int dim>
+double RightHandSide<dim>::value (const Point<dim>  &/*p*/,
+                                 const unsigned int /*component*/) const 
+{
+  return 0;
+}
+
+
+
+template <int dim>
+double PressureBoundaryValues<dim>::value (const Point<dim>  &p,
+                                          const unsigned int /*component*/) const 
+{
+  return 1-p[0];
+}
+
+
+template <int dim>
+void
+SaturationBoundaryValues<dim>::vector_value (const Point<dim> &p,
+                                            Vector<double>   &values) const 
+{
+  Assert (values.size() == dim+2,
+         ExcDimensionMismatch (values.size(), dim+2));
+
+  for (unsigned int d=0; d<dim+1; ++d)
+    values(d) = 0;
+
+  if (p[0] == 0)
+    values(dim+1) = 1;
+  else
+    values(dim+1) = 0;
+}
+
+
+
+template <int dim>
+double InitialValues<dim>::value (const Point<dim>  &p,
+                                          const unsigned int component) const 
+{
+    if(component<dim+1)
+    return 0;
+    else 
+    { 
+      if(p[0]==0)return 1;
+      else return 0;
+    }
+  
+}
+
+
+template <int dim>
+void
+InitialValues<dim>::vector_value (const Point<dim> &p,
+                                 Vector<double>   &values) const 
+{
+  Assert (values.size() == dim+2,
+         ExcDimensionMismatch (values.size(), dim+2));
+
+  for (unsigned int d=0; d<dim+1; ++d)
+    values(d) = 0;
+  values(dim+1) = InitialValues::value(p,dim+1);
+}
+
+
+
+
+                                // {The inverse permeability tensor, 
+                                      Coefficient and inverse mobility scalar}
+
+                                
+                                 //For the inverse  permeability tensor,
+                                 // ``KInverse''.As in introduction, '
+                                 // assume the heterogeneous is isotropic,
+                                 // so it is a scalar multipy the identity matrix.
+                                  //DealII has a base class not only for
+                                 // scalar and generally vector-valued
+                                 // functions (the ``Function'' base
+                                 // class) but also for functions that
+                                 // return tensors of fixed dimension
+                                 // and rank, the ``TensorFunction''
+                                 // template. Here, the function under
+                                 // consideration returns a dim-by-dim
+                                 // matrix, i.e. a tensor of rank 2
+                                 // and dimension ``dim''. We then
+                                 // choose the template arguments of
+                                 // the base class appropriately.
+                                 //
+                                 // The interface that the
+                                 // ``TensorFunction'' class provides
+                                 // is essentially equivalent to the
+                                 // ``Function'' class. In particular,
+                                 // there exists a ``value_list''
+                                 // function that takes a list of
+                                 // points at which to evaluate the
+                                 // function, and returns the values
+                                 // of the function in the second
+                                 // argument, a list of tensors:
+template <int dim>
+class KInverse : public TensorFunction<2,dim>
+{
+  public:
+    virtual void value_list (const std::vector<Point<dim> > &points,
+                            std::vector<Tensor<2,dim> >    &values,
+                             const double epsilon) const;
+};
+
+
+template <int dim>
+class Coefficient : public Function<dim>
+{
+  public:
+    Coefficient () : Function<dim>(1) {};
+    
+    virtual double value (const Point<dim>   &p,
+                         const unsigned int  component = 0 ) const;
+};
+
+
+template <int dim>
+void
+KInverse<dim>::value_list (const std::vector<Point<dim> > &points,
+                           std::vector<Tensor<2,dim> >    &values,
+                           const double epsilon) const
+{
+  Assert (points.size() == values.size(),
+         ExcDimensionMismatch (points.size(), values.size()));
+  for (unsigned int p=0; p<points.size(); ++p)
+    {
+      values[p].clear ();
+
+      const double permeability = 1.0/(2+1.99*sin(2*3.1415926*(2*points[p][0]-1*points[p][1])/epsilon));
+      
+      for (unsigned int d=0; d<dim; ++d)
+          values[p][d][d]=1.0/permeability;
+
+    }
+}
+
+
+double mobility_inverse (const double S, const double vis)
+{ 
+   return 1.0 /(1.0/vis * S * S + (1-S) * (1-S));
+}
+
+double f_saturation(const double S, const double vis)
+{   
+
+   return S*S /( S * S +vis * (1-S) * (1-S));
+}
+
+template <int dim>
+double Coefficient<dim>::value (const Point<dim>  &/*p*/,
+                                 const unsigned int /*component*/) const 
+{
+  return 1;
+}
+
+
+
+
+                                 // {extract_u and friends}
+
+                                 // The next five functions are
+                                 // needed for matrix and right hand
+                                 // side assembly. They are described
+                                 // in detail in step-20:
+template <int dim>
+Tensor<1,dim>
+extract_u (const FEValuesBase<dim> &fe_values,
+          const unsigned int i,
+          const unsigned int q)
+{
+  Tensor<1,dim> tmp;
+
+  for (unsigned int d=0; d<dim; ++d)
+    tmp[d] = fe_values.shape_value_component (i,q,d);
+
+  return tmp;
+}
+
+
+
+template <int dim>
+double
+extract_div_u (const FEValuesBase<dim> &fe_values,
+              const unsigned int i,
+              const unsigned int q)
+{
+  double divergence = 0;
+  for (unsigned int d=0; d<dim; ++d)
+    divergence += fe_values.shape_grad_component (i,q,d)[d];
+
+  return divergence;
+}
+
+
+  
+template <int dim>
+double extract_p (const FEValuesBase<dim> &fe_values,
+                  const unsigned int i,
+                  const unsigned int q)
+{
+  return fe_values.shape_value_component (i,q,dim);
+}
+
+template <int dim>
+double extract_s (const FEValuesBase<dim> &fe_values,
+                  const unsigned int i,
+                  const unsigned int q)
+{
+  return fe_values.shape_value_component (i,q,dim+1);
+}
+
+template <int dim>
+Tensor<1,dim>
+extract_grad_s(const FEValuesBase<dim> &fe_values,
+              const unsigned int i,
+              const unsigned int q)
+{
+  Tensor<1,dim> tmp;
+  for (unsigned int d=0; d<dim; ++d)
+    tmp[d] = fe_values.shape_grad_component (i,q,dim+1)[d];
+
+  return tmp;
+}
+
+
+
+                                 // {TwoPhaseFlowProblem class implementation}
+
+                                 // {TwoPhaseFlowProblem::TwoPhaseFlowProblem}
+                                 //  we use RT(k) X DG(k),DG(k) spaces.
+                                 // time_step is small enough to make the solution 
+                                 // converges stably. 
+
+                               
+                                 
+template <int dim>
+TwoPhaseFlowProblem<dim>::TwoPhaseFlowProblem (const unsigned int degree)
+               :
+               degree (degree),
+                fe (FE_RaviartThomas<dim>(degree), 1,
+                    FE_DGQ<dim>(degree), 1,
+                   FE_DGQ<dim>(degree), 1),
+               dof_handler (triangulation),
+               n_refinement_steps (5),
+               time_step (1.0/std::pow(2.0, double(n_refinement_steps))/6),
+                epsilon(0.05),
+                vis (0.2)
+                
+{}
+
+
+
+                                 // {TwoPhaseFlowProblem::make_grid_and_dofs}
+
+                                 // This next function starts out with
+                                 // well-known functions calls that
+                                 // create and refine a mesh, and then
+                                 // associate degrees of freedom with
+                                 // it:
+template <int dim>
+void TwoPhaseFlowProblem<dim>::make_grid_and_dofs ()
+{
+  GridGenerator::hyper_cube (triangulation, 0, 1);
+  
+  for (unsigned int f=0; f<GeometryInfo<dim>::faces_per_cell; ++f)
+   { if (triangulation.begin()->face(f)->center()[0] == 0)
+      triangulation.begin()->face(f)->set_boundary_indicator (1);
+     if (triangulation.begin()->face(f)->center()[0] == 1)
+      triangulation.begin()->face(f)->set_boundary_indicator (2);
+   }
+
+  triangulation.refine_global (n_refinement_steps);
+  
+  dof_handler.distribute_dofs (fe);
+
+  DoFRenumbering::component_wise (dof_handler);
+
+                                  
+  std::vector<unsigned int> dofs_per_component (dim+2);
+  DoFTools::count_dofs_per_component (dof_handler, dofs_per_component);  
+  const unsigned int n_u = dofs_per_component[0],
+                     n_p = dofs_per_component[dim],
+                    n_s = dofs_per_component[dim+1];
+
+  std::cout << "Number of active cells: "
+           << triangulation.n_active_cells()
+           << std::endl
+           << "Total number of cells: "
+           << triangulation.n_cells()
+           << std::endl
+            << "Number of degrees of freedom: "
+           << dof_handler.n_dofs()
+            << " (" << n_u << '+' << n_p << '+'<< n_s <<')'
+           << std::endl;
+
+  
+  const unsigned int
+    n_couplings = dof_handler.max_couplings_between_dofs();
+  
+  sparsity_pattern.reinit (3,3);
+  sparsity_pattern.block(0,0).reinit (n_u, n_u, n_couplings);
+  sparsity_pattern.block(1,0).reinit (n_p, n_u, n_couplings);
+  sparsity_pattern.block(2,0).reinit (n_s, n_u, n_couplings);
+  sparsity_pattern.block(0,1).reinit (n_u, n_p, n_couplings);
+  sparsity_pattern.block(1,1).reinit (n_p, n_p, n_couplings);
+  sparsity_pattern.block(2,1).reinit (n_s, n_p, n_couplings);
+  sparsity_pattern.block(0,2).reinit (n_u, n_s, n_couplings);
+  sparsity_pattern.block(1,2).reinit (n_p, n_s, n_couplings);
+  sparsity_pattern.block(2,2).reinit (n_s, n_s, n_couplings);
+  
+  sparsity_pattern.collect_sizes();
+
+  
+  DoFTools::make_sparsity_pattern (dof_handler, sparsity_pattern);
+  sparsity_pattern.compress();
+
+  system_matrix.reinit (sparsity_pattern);
+
+                                   
+  solution.reinit (3);
+  solution.block(0).reinit (n_u);
+  solution.block(1).reinit (n_p);
+  solution.block(2).reinit (n_s);
+  solution.collect_sizes ();
+  
+  old_solution.reinit (3);
+  old_solution.block(0).reinit (n_u);
+  old_solution.block(1).reinit (n_p);
+  old_solution.block(2).reinit (n_s);
+  old_solution.collect_sizes ();
+  
+  system_rhs.reinit (3);
+  system_rhs.block(0).reinit (n_u);
+  system_rhs.block(1).reinit (n_p);
+  system_rhs.block(2).reinit (n_s);
+  system_rhs.collect_sizes ();
+
+
+}
+
+
+                                 // {TwoPhaseFlowProblem::assemble_system}
+                                 // The function that
+                                 // assembles the linear system has
+                                 // mostly been discussed already in
+                                 // the introduction to this
+                                 // test case. We want to emphasize that
+                                 // we assemble the first two equations
+                                 // for velocity and pressure, but 
+                                 // for saturation we only assemble 
+                                 // the Matrixblock(2,2), for Matrixblock(0,2)
+                                 // we will assemble it in "solve()", because
+                                 //at that time, we have the new velocity solved
+                                 // we can use it to assemble Matrixblock(0,2)
+                    
+
+template <int dim>
+void TwoPhaseFlowProblem<dim>::assemble_system () 
+{  
+  QGauss<dim>   quadrature_formula(degree+2); 
+  QGauss<dim-1> face_quadrature_formula(degree+2);
+
+  FEValues<dim> fe_values (fe, quadrature_formula, 
+                          update_values    | update_gradients |
+                           update_q_points  | update_JxW_values);
+  FEFaceValues<dim> fe_face_values (fe, face_quadrature_formula, 
+                                   update_values    | update_normal_vectors |
+                                   update_q_points  | update_JxW_values);
+
+  const unsigned int   dofs_per_cell   = fe.dofs_per_cell;
+  
+  const unsigned int   n_q_points      = quadrature_formula.n_quadrature_points;
+  const unsigned int   n_face_q_points = face_quadrature_formula.n_quadrature_points;
+
+  FullMatrix<double>   local_matrix (dofs_per_cell, dofs_per_cell);
+  Vector<double>       local_rhs (dofs_per_cell);
+
+  std::vector<unsigned int> local_dof_indices (dofs_per_cell);
+  
+                                   // The next step is to declare
+                                   // objects that represent the
+                                   // source term, pressure boundary
+                                   // value, and coefficient in the
+                                   // equation. In addition to these
+                                   // objects that represent
+                                   // continuous functions, we also
+                                   // need arrays to hold their values
+                                   // at the quadrature points of
+                                   // individual cells (or faces, for
+                                   // the boundary values). Note that
+                                   // in the case of the coefficient,
+                                   // the array has to be one of
+                                   // matrices.
+  const RightHandSide<dim>          right_hand_side;
+  const PressureBoundaryValues<dim> pressure_boundary_values;
+  const KInverse<dim>               k_inverse;
+  const Coefficient<dim>            coefficient;
+   
+  
+  std::vector<double>               rhs_values (n_q_points);
+  std::vector<double>               boundary_values (n_face_q_points);
+  std::vector<Tensor<2,dim> >       k_inverse_values (n_q_points);
+  std::vector<double>               coefficient_values(n_q_points);
+  
+  std::vector<Vector<double> >      old_solution_values(n_q_points, Vector<double>(dim+2));
+  std::vector<std::vector<Tensor<1,dim> > >  old_solution_grads(n_q_points,
+                                                                std::vector<Tensor<1,dim> > (dim+2));
+  
+  
+
+                                   // With all this in place, we can
+                                   // go on with the loop over all
+                                   // cells. The body of this loop has
+                                   // been discussed in the
+                                   // introduction, and will not be
+                                   // commented any further here:
+  typename DoFHandler<dim>::active_cell_iterator
+    cell = dof_handler.begin_active(),
+    endc = dof_handler.end();
+  unsigned int cellnum=0;
+  system_matrix=0;
+  system_rhs=0;
+  for (; cell!=endc; ++cell)
+    { cellnum++;
+      fe_values.reinit (cell);
+      local_matrix = 0;
+      local_rhs = 0;
+
+      fe_values.get_function_values (old_solution, old_solution_values);
+      right_hand_side.value_list (fe_values.get_quadrature_points(),
+                                  rhs_values);
+      k_inverse.value_list (fe_values.get_quadrature_points(),
+                            k_inverse_values,
+                            epsilon);
+      
+      coefficient.value_list (fe_values.get_quadrature_points(), coefficient_values);
+      
+      for (unsigned int q=0; q<n_q_points; ++q)            
+        for (unsigned int i=0; i<dofs_per_cell; ++i)
+          {
+           const double old_s = old_solution_values[q](dim+1);
+
+            const Tensor<1,dim> phi_i_u = extract_u (fe_values, i, q);
+           const double div_phi_i_u = extract_div_u (fe_values, i, q);
+            const double phi_i_p = extract_p (fe_values, i, q);
+           const double phi_i_s = extract_s (fe_values, i, q); 
+           const Tensor<1,dim> grad_phi_i_s = extract_grad_s(fe_values, i, q);
+                    
+            
+            for (unsigned int j=0; j<dofs_per_cell; ++j)
+              {
+                const Tensor<1,dim> phi_j_u = extract_u (fe_values, j, q);
+               const double div_phi_j_u = extract_div_u (fe_values, j, q);             
+                const double phi_j_p = extract_p (fe_values, j, q);
+                const double phi_j_s = extract_s (fe_values, j, q);  
+               
+                local_matrix(i,j) += (phi_i_u * k_inverse_values[q] *
+                                     mobility_inverse(old_s,vis) * phi_j_u            
+                                      - div_phi_i_u * phi_j_p
+                                      - phi_i_p * div_phi_j_u
+                                     + coefficient_values[q] * phi_i_s * phi_j_s
+                                     )
+                                     * fe_values.JxW(q);     
+              }
+
+            local_rhs(i) += (-phi_i_p * rhs_values[q])*
+                            fe_values.JxW(q);
+          }
+      
+                                             //here, we compute the boundary values for pressure 
+
+      for (unsigned int face_no=0;
+          face_no<GeometryInfo<dim>::faces_per_cell;
+          ++face_no)
+       if (cell->at_boundary(face_no))
+         {
+           fe_face_values.reinit (cell, face_no);
+           
+           pressure_boundary_values
+             .value_list (fe_face_values.get_quadrature_points(),
+                          boundary_values);
+
+           for (unsigned int q=0; q<n_face_q_points; ++q) 
+             for (unsigned int i=0; i<dofs_per_cell; ++i)
+               {
+                 const Tensor<1,dim>
+                   phi_i_u = extract_u (fe_face_values, i, q);
+
+                 local_rhs(i) += -(phi_i_u *
+                                   fe_face_values.normal_vector(q) *
+                                   boundary_values[q] *
+                                   fe_face_values.JxW(q));
+               }
+         }
+
+                                       // The final step in the loop
+                                       // over all cells is to
+                                       // transfer local contributions
+                                       // into the global matrix and
+                                       // right hand side vector. Note
+                                       // that we use exactly the same
+                                       // interface as in previous
+                                       // examples, although we now
+                                       // use block matrices and
+                                       // vectors instead of the
+                                       // regular ones. In other
+                                       // words, to the outside world,
+                                       // block objects have the same
+                                       // interface as matrices and
+                                       // vectors, but they
+                                       // additionally allow to access
+                                       // individual blocks.
+      cell->get_dof_indices (local_dof_indices);
+      for (unsigned int i=0; i<dofs_per_cell; ++i)
+      
+        for (unsigned int j=0; j<dofs_per_cell; ++j)
+      {    system_matrix.add (local_dof_indices[i],
+                             local_dof_indices[j],
+                             local_matrix(i,j));
+       }
+      
+      for (unsigned int i=0; i<dofs_per_cell; ++i)
+        system_rhs(local_dof_indices[i]) += local_rhs(i);      
+       
+    }
+}
+
+
+                                 // {Linear solvers and preconditioners}
+
+                                 // {The ``InverseMatrix'' class template}
+                                 
+                                // Everything here is completely same with step-20
+                                 
+
+
+template <class Matrix>
+class InverseMatrix : public Subscriptor
+{
+  public:
+    InverseMatrix (const Matrix &m);
+
+    void vmult (Vector<double>       &dst,
+                const Vector<double> &src) const;
+
+  private:
+    const SmartPointer<const Matrix> matrix;
+
+    mutable GrowingVectorMemory<> vector_memory;    
+};
+
+
+template <class Matrix>
+InverseMatrix<Matrix>::InverseMatrix (const Matrix &m)
+                :
+                matrix (&m)
+{}
+
+
+                                 
+template <class Matrix>
+void InverseMatrix<Matrix>::vmult (Vector<double>       &dst,
+                                   const Vector<double> &src) const
+{
+  SolverControl solver_control (src.size(), 1e-8*src.l2_norm());
+  SolverCG<> cg (solver_control, vector_memory);
+
+  dst = 0;
+  
+  cg.solve (*matrix, dst, src, PreconditionIdentity());        
+}
+
+
+                                 // {The ``SchurComplement'' class template}
+                                                                 
+class SchurComplement : public Subscriptor
+{
+  public:
+    SchurComplement (const BlockSparseMatrix<double> &A,
+                     const InverseMatrix<SparseMatrix<double> > &Minv);
+
+    void vmult (Vector<double>       &dst,
+                const Vector<double> &src) const;
+
+  private:
+    const SmartPointer<const BlockSparseMatrix<double> > system_matrix;
+    const SmartPointer<const InverseMatrix<SparseMatrix<double> > > m_inverse;
+    
+    mutable Vector<double> tmp1, tmp2;
+};
+
+
+SchurComplement::SchurComplement (const BlockSparseMatrix<double> &A,
+                                  const InverseMatrix<SparseMatrix<double> > &Minv)
+                :
+                system_matrix (&A),
+                m_inverse (&Minv),
+                tmp1 (A.block(0,0).m()),
+                tmp2 (A.block(0,0).m())
+{}
+
+
+void SchurComplement::vmult (Vector<double>       &dst,
+                             const Vector<double> &src) const
+{
+  system_matrix->block(0,1).vmult (tmp1, src);
+  m_inverse->vmult (tmp2, tmp1);
+  system_matrix->block(1,0).vmult (dst, tmp2);
+}
+
+
+                                 // {The ``ApproximateSchurComplement'' class template}
+
+class ApproximateSchurComplement : public Subscriptor
+{
+  public:
+    ApproximateSchurComplement (const BlockSparseMatrix<double> &A);
+
+    void vmult (Vector<double>       &dst,
+                const Vector<double> &src) const;
+
+  private:
+    const SmartPointer<const BlockSparseMatrix<double> > system_matrix;
+    
+    mutable Vector<double> tmp1, tmp2;
+};
+
+
+ApproximateSchurComplement::ApproximateSchurComplement (const BlockSparseMatrix<double> &A)
+                :
+                system_matrix (&A),
+                tmp1 (A.block(0,0).m()),
+                tmp2 (A.block(0,0).m())
+{}
+
+
+void ApproximateSchurComplement::vmult (Vector<double>       &dst,
+                                        const Vector<double> &src) const
+{
+  system_matrix->block(0,1).vmult (tmp1, src);
+  system_matrix->block(0,0).precondition_Jacobi (tmp2, tmp1);
+  system_matrix->block(1,0).vmult (dst, tmp2);
+}
+
+
+
+                                 // {TwoPhaseFlowProblem::solve}
+
+                                 // After all these preparations,
+                                 // we finally solves the linear
+                                 // system for velocity and pressure.
+                                 // And remember, we still have to assemble 
+                                 // the Matirxbloc(2,0) after velocity is computed
+                                 // , then use it to solve saturation.
+template <int dim>
+void TwoPhaseFlowProblem<dim>::solve () 
+{
+  const InverseMatrix<SparseMatrix<double> >
+    m_inverse (system_matrix.block(0,0));
+  Vector<double> tmp (solution.block(0).size());
+  Vector<double> schur_rhs (solution.block(1).size());
+  Vector<double> tmp2 (solution.block(2).size());
+  
+
+                                 // this part is for pressure
+  {
+    m_inverse.vmult (tmp, system_rhs.block(0));
+    system_matrix.block(1,0).vmult (schur_rhs, tmp);
+    schur_rhs -= system_rhs.block(1);
+
+    
+    SchurComplement
+      schur_complement (system_matrix, m_inverse);
+    
+    ApproximateSchurComplement
+      approximate_schur_complement (system_matrix);
+      
+    InverseMatrix<ApproximateSchurComplement>
+      preconditioner (approximate_schur_complement);
+
+    
+    SolverControl solver_control (system_matrix.block(0,0).m(),
+                                 1e-12*schur_rhs.l2_norm());
+    SolverCG<>    cg (solver_control);
+
+    cg.solve (schur_complement, solution.block(1), schur_rhs,
+              preconditioner);
+  
+    std::cout << solver_control.last_step()
+              << " CG Schur complement iterations to obtain convergence for pressure."
+              << std::endl;
+  }
+
+                                   //  this part is for velocity. The
+                                   // equation reads MU=-B^TP+F, and
+                                   // we solve it by first computing
+                                   // the right hand side, and then
+                                   // multiplying it with the object
+                                   // that represents the inverse of
+                                   // the mass matrix:
+  {
+    system_matrix.block(0,1).vmult (tmp, solution.block(1));
+    tmp *= -1;
+    tmp += system_rhs.block(0);
+    
+    m_inverse.vmult (solution.block(0), tmp);
+  }
+
+                                      //This part is for saturation.
+                                      // Here are many complicated functions
+                                      //which are very similiar with the
+                                      //assemble_system() part.
+                                     // For DG(0), we have to consider the discontinuty
+                                     // of the solution, then as in Introduction,
+                                     // compute numerical flux and judge it is in-flow or out-flow.
+                                     // After assemble Matrixbloc(2,0)
+                                     // , we could compute saturation directly. 
+  { 
+  QGauss<dim>   quadrature_formula(degree+2); 
+  QGauss<dim-1> face_quadrature_formula(degree+2);  
+  FEValues<dim> fe_values (fe, quadrature_formula, 
+                          update_values    | update_gradients |
+                           update_q_points  | update_JxW_values);
+  FEFaceValues<dim> fe_face_values (fe, face_quadrature_formula, 
+                                   update_values    | update_normal_vectors |
+                                   update_q_points  | update_JxW_values);
+  FEFaceValues<dim> fe_face_values_neighbor (fe, face_quadrature_formula, 
+                                            update_values);
+  
+  const unsigned int   dofs_per_cell   = fe.dofs_per_cell;
+  const unsigned int   n_q_points      = quadrature_formula.n_quadrature_points;
+  const unsigned int   n_face_q_points = face_quadrature_formula.n_quadrature_points;
+  
+  vfs_out = 0.0;
+  v_out = 0.0;  
+  
+  Vector<double>       local_rhs (dofs_per_cell);
+  std::vector<Vector<double> > old_solution_values(n_q_points, Vector<double>(dim+2));
+  std::vector<Vector<double> > old_solution_values_face(n_face_q_points, Vector<double>(dim+2));
+  std::vector<Vector<double> > old_solution_values_face_neighbor(n_face_q_points, Vector<double>(dim+2));
+  std::vector<Vector<double> > present_solution_values(n_q_points, Vector<double>(dim+2));
+  std::vector<Vector<double> > present_solution_values_face(n_face_q_points, Vector<double>(dim+2));
+
+  std::vector<double> neighbor_saturation (n_face_q_points);
+  std::vector<unsigned int> local_dof_indices (dofs_per_cell);
+  
+  typename DoFHandler<dim>::active_cell_iterator
+    cell = dof_handler.begin_active(),
+    endc = dof_handler.end();
+
+ for (; cell!=endc; ++cell)
+ {
+    local_rhs = 0;
+    fe_values.reinit (cell);
+
+    fe_values.get_function_values (old_solution, old_solution_values);
+    fe_values.get_function_values (solution, present_solution_values);
+    
+    for (unsigned int q=0; q<n_q_points; ++q) 
+      for (unsigned int i=0; i<dofs_per_cell; ++i)
+       {
+         const double old_s = old_solution_values[q](dim+1);
+          Tensor<1,dim> present_u;
+         for (unsigned int d=0; d<dim; ++d)
+             present_u[d] = present_solution_values[q](d);
+
+         const double phi_i_s = extract_s(fe_values, i, q);
+         const Tensor<1,dim> grad_phi_i_s = extract_grad_s(fe_values, i, q);
+                    
+         local_rhs(i) += (
+                           time_step *(f_saturation(old_s,vis) * present_u * grad_phi_i_s)+
+                            old_s * phi_i_s)
+                         * fe_values.JxW(q);
+               }
+                                       //Here is our numerical flux computation
+                                        // Finding neighbor as step-12
+                                 
+    for (unsigned int face_no=0; face_no<GeometryInfo<dim>::faces_per_cell;++face_no)
+       {
+        fe_face_values.reinit (cell, face_no);
+
+        fe_face_values.get_function_values (old_solution, old_solution_values_face);
+        fe_face_values.get_function_values (solution, present_solution_values_face);
+
+        if (cell->at_boundary(face_no))
+          {
+            if (cell->face(face_no)->boundary_indicator() == 1)
+              for (unsigned int q=0;q<n_face_q_points;++q)
+                neighbor_saturation[q] = 1;
+            else
+              for (unsigned int q=0;q<n_face_q_points;++q)
+                neighbor_saturation[q] = 0;                     
+          }
+        else
+                                           // there is a neighbor behind this face
+          {
+            const typename DoFHandler<dim>::active_cell_iterator
+              neighbor = cell->neighbor(face_no);
+            const unsigned int
+              neighbor_face = cell->neighbor_of_neighbor(face_no);
+
+            fe_face_values_neighbor.reinit (neighbor, neighbor_face);
+            
+            fe_face_values_neighbor.get_function_values (old_solution,
+                                                         old_solution_values_face_neighbor);
+            
+            for (unsigned int q=0;q<n_face_q_points;++q)
+              neighbor_saturation[q] = old_solution_values_face_neighbor[q](dim+1);
+          }
+          
+
+         if (cell->at_boundary(face_no))
+           {   
+              if (cell->face(face_no)->boundary_indicator() ==2 )
+               {for (unsigned int q=0;q<n_face_q_points;++q)
+                 {
+                    vfs_out += present_solution_values_face[q](0)
+                             *f_saturation(present_solution_values_face[q](dim+1),vis)
+                             *fe_face_values.JxW(q);
+                     v_out += present_solution_values_face[q](0)
+                             *fe_face_values.JxW(q);
+                   }                                
+                }
+             }
+        for (unsigned int q=0;q<n_face_q_points;++q)
+          {
+            Tensor<1,dim> present_u_face;
+            for (unsigned int d=0; d<dim; ++d)
+             { present_u_face[d] = present_solution_values_face[q](d);
+     }
+            const double normal_flux = present_u_face *
+                                       fe_face_values.normal_vector(q);
+
+            const bool is_outflow_q_point = (normal_flux >= 0);
+                                    
+            if (is_outflow_q_point == true)
+              {
+                for (unsigned int i=0; i<dofs_per_cell; ++i)
+                { 
+                   const double outflow = -time_step * normal_flux 
+                                        * f_saturation(old_solution_values_face[q](dim+1),vis)
+                                        * extract_s(fe_face_values,i,q)
+                                        * fe_face_values.JxW(q);
+                   local_rhs(i) += outflow;
+                } 
+              }
+             
+            else
+             {
+               for (unsigned int i=0; i<dofs_per_cell; ++i)
+                {
+                  const double inflow = -time_step * normal_flux 
+                                       * f_saturation( neighbor_saturation[q],vis)
+                                       * extract_s(fe_face_values,i,q)
+                                        * fe_face_values.JxW(q);
+                 local_rhs(i) += inflow;
+                }
+               
+            }
+       
+          }
+             
+       }
+  
+     cell->get_dof_indices (local_dof_indices);
+     for (unsigned int i=0; i<dofs_per_cell; ++i)
+       {
+         system_rhs(local_dof_indices[i]) += local_rhs(i);
+        }
+               
+   }   
+    SolverControl solver_control (system_matrix.block(2,2).m(),
+                                 1e-12*system_rhs.block(2).l2_norm());
+    SolverCG<>   cg (solver_control);
+    cg.solve (system_matrix.block(2,2), solution.block(2), system_rhs.block(2),
+               PreconditionIdentity());
+               
+       
+    std::cout << solver_control.last_step()
+              << " CG iterations to obtain convergence for saturation."
+              << std::endl;            
+  } 
+
+   
+    old_solution = solution; 
+
+   
+  
+}
+                                 
+                                 // {TwoPhaseFlow::compute_errors}
+
+                                 // After we have dealt with the
+                                 // linear solver and preconditioners,
+                                 // we continue with the
+                                 // implementation of our main
+                                 // class. In particular, the next
+                                 // task is to compute the errors in
+                                 // our numerical solution, in both
+                                 // the pressures velocities as well as
+                                 // saturations.
+                                 //
+                                 // To compute errors in the solution,
+                                 // we will not use ``VectorTools::integrate_difference''
+                                 // as  step-20,  since we don't have exact solutions.
+                                 // What we will do is to give some points
+                                 // and evaluate the values on these points.
+                                 //For every solution, we get values on those points,
+                                 // then we can compare the values as an error.
+   
+
+template <int dim>
+Vector<double>
+TwoPhaseFlowProblem<dim>::evaluate_solution (const Point<dim> &point) const
+{
+  static const MappingQ1<dim> mapping;
+                                    // first find the cell in which this point
+                                            // is, initialize a quadrature rule with
+                                   // it, and then a FEValues object
+  const typename DoFHandler<dim>::active_cell_iterator
+    cell = GridTools::find_active_cell_around_point (dof_handler, point);
+
+  const Point<dim> unit_point
+    = mapping.transform_real_to_unit_cell(cell, point);
+  Assert (GeometryInfo<dim>::is_inside_unit_cell (unit_point),
+          ExcInternalError());
+
+  const Quadrature<dim> quadrature (unit_point);
+  FEValues<dim> fe_values(mapping, fe, quadrature, update_values);
+  fe_values.reinit(cell);
+                                   // then use this to get at the values of
+                                   // the given fe_function at this point
+  std::vector<Vector<double> > u_value(1, Vector<double>(dim+2));
+  fe_values.get_function_values(solution, u_value);
+
+  return u_value[0];
+}
+
+                                  //{TwoPhaseFlowProblem::compute_errors}
+
+                                  // The compute_errors function is to compute
+                                  // error on some euqally spaced fixed points
+                                  // use evaluation function to interpret 
+                                  // solution value at the point
+                                  // then output those fixed points' value
+                                  // For each mesh, we can compare the output
+                                  // to estimate errors.
+   
+template <int dim>
+void TwoPhaseFlowProblem<dim>::compute_errors () const
+{
+  std::ofstream sampled_solution ("sampled_solution");
+
+  const double dx = 0.01;
+  const double dy = 0.01;
+
+  for (double x=0; x<=1; x+=dx)
+    for (double y=0; y<=1; y+=dy)
+      {
+       const Point<dim> point(x,y);
+
+       Vector<double> solution_at_point(dim+2);
+
+       solution_at_point = evaluate_solution (point);
+
+       sampled_solution << point << " ";
+       for (unsigned int c=0; c<dim+2; ++c)
+         sampled_solution << solution_at_point(c) << " ";
+       sampled_solution << std::endl;
+      }
+}
+
+                                 // {TwoPhaseFlowProblem::output_results}
+
+                                 // The output_results function is
+                                 // the one in which we generate
+                                 // graphical output.
+template <int dim>
+void TwoPhaseFlowProblem<dim>::output_results 
+(const unsigned int timestep_number)  const
+{  
+  std::vector<std::string> solution_names;
+  switch (dim)
+    {
+      case 2:
+            solution_names.push_back ("u");
+            solution_names.push_back ("v");
+            solution_names.push_back ("p");
+           solution_names.push_back ("S");
+            break;
+            
+      case 3:
+            solution_names.push_back ("u");
+            solution_names.push_back ("v");
+            solution_names.push_back ("w");
+            solution_names.push_back ("p");
+           solution_names.push_back ("S");
+            break;
+            
+      default:
+            Assert (false, ExcNotImplemented());
+    }
+  
+  DataOut<dim> data_out;
+
+  data_out.attach_dof_handler (dof_handler);
+  data_out.add_data_vector (solution, solution_names);
+
+  data_out.build_patches (degree+1);
+  
+  std::ostringstream filename;
+  filename << "solution-"<< timestep_number;
+
+  std::ofstream output (filename.str().c_str());
+  data_out.write_gnuplot (output);
+
+  //data_out.write_vtk (output);
+}
+
+
+                                 // {TwoPhaseFlowProblem::run}
+
+                                 // This is the final function of our
+                                 // main class. It's only job is to
+                                 // call the other functions in their order:
+template <int dim>
+void TwoPhaseFlowProblem<dim>::run () 
+{
+  std::cout<<"Solving problem in " <<dim << " space dimensions." << std::endl;
+  
+  make_grid_and_dofs();
+  
+  ConstraintMatrix constraints;
+  constraints.close();
+
+  std::list<double> production_rate;
+  std::list<double> production_time;
+
+  Vector<double> tmp (old_solution.size());
+  VectorTools::project (dof_handler, constraints, QGauss<dim>(degree+2),InitialValues<dim>(),tmp);
+  std::copy (tmp.begin(), tmp.end(), old_solution.begin());
+  
+  unsigned int timestep_number = 1;
+  
+  for ( double time = time_step; time <=1; time+=time_step,  timestep_number++)
+  { 
+    std::cout<< "Timestep_number = "<< timestep_number<<std::endl; 
+    assemble_system ();
+    solve ();
+    output_results(timestep_number);
+
+    production_time.push_back (time);
+    production_rate.push_back (1.0 - vfs_out/v_out);
+    std::cout<<"production_rate="<<production_rate.back()<<std::endl;       
+  }
+
+  std::ofstream production_history ("production_history");
+  std::list<double>::iterator
+    list_element = production_rate.begin(),
+    time_element = production_time.begin();
+  for (; list_element != production_rate.end(); ++list_element, ++time_element)
+    production_history << *time_element << " " << *list_element << std::endl;
+  
+   
+  compute_errors ();
+}
+
+    
+                                 // {The ``main'' function}
+
+                                // In the main function, we pass the
+                                // degree of the finite element space
+                                // to the constructor of the TwoPhaseFlowProblem
+                                // (here, we use zero-th order elements).
+int main () 
+{
+  try
+    {
+      deallog.depth_console (0);
+
+      TwoPhaseFlowProblem<2> two_phase_flow_problem(0);
+      two_phase_flow_problem.run ();
+    }
+  catch (std::exception &exc)
+    {
+      std::cerr << std::endl << std::endl
+               << "----------------------------------------------------"
+               << std::endl;
+      std::cerr << "Exception on processing: " << std::endl
+               << exc.what() << std::endl
+               << "Aborting!" << std::endl
+               << "----------------------------------------------------"
+               << std::endl;
+      
+      return 1;
+    }
+  catch (...) 
+    {
+      std::cerr << std::endl << std::endl
+               << "----------------------------------------------------"
+               << std::endl;
+      std::cerr << "Unknown exception!" << std::endl
+               << "Aborting!" << std::endl
+               << "----------------------------------------------------"
+               << std::endl;
+      return 1;
+    }
+
+  return 0;
+}

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.