From: bangerth Date: Tue, 17 Oct 2006 04:24:27 +0000 (+0000) Subject: Add Yan Li's step-21 for the IMPES problem X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=b9cea7e1c41e35521389d5bf00c70e535fe43b40;p=dealii-svn.git Add Yan Li's step-21 for the IMPES problem git-svn-id: https://svn.dealii.org/trunk@14010 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/examples/step-21/Makefile b/deal.II/examples/step-21/Makefile new file mode 100644 index 0000000000..1972f5463f --- /dev/null +++ b/deal.II/examples/step-21/Makefile @@ -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========= $( $@ \ + || (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 index 0000000000..f693a95d3c --- /dev/null +++ b/deal.II/examples/step-21/doc/intro.dox @@ -0,0 +1 @@ +

Introduction

diff --git a/deal.II/examples/step-21/doc/perm.png b/deal.II/examples/step-21/doc/perm.png new file mode 100644 index 0000000000..1c7fe71ebc 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 index 0000000000..0f0afe5798 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 index 0000000000..89faa40ac4 --- /dev/null +++ b/deal.II/examples/step-21/doc/results.dox @@ -0,0 +1 @@ +

Results

diff --git a/deal.II/examples/step-21/doc/solp4meshes.png b/deal.II/examples/step-21/doc/solp4meshes.png new file mode 100644 index 0000000000..1a15adf7c3 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 index 0000000000..9cce4b42cd 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 index 0000000000..8264339c3c 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 index 0000000000..4b1971786c 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 index 0000000000..64779eb700 --- /dev/null +++ b/deal.II/examples/step-21/doc/yanli_project_introduction.tex @@ -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 index 0000000000..699f6e6990 --- /dev/null +++ b/deal.II/examples/step-21/doc/yanli_project_result.tex @@ -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 index 0000000000..ecb30a7416 --- /dev/null +++ b/deal.II/examples/step-21/step-21.cc @@ -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 +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + //The Discontinuous Galerkin finite element is declared: +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include + // The Raviart-Thomas finite element is declared: +#include + + // 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 + + + // {The ``TwoPhaseFlowProblem'' class template} + + +template +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 evaluate_solution (const Point &point) const; + + const unsigned int degree; + + Triangulation triangulation; + FESystem fe; + DoFHandler dof_handler; + + + BlockSparsityPattern sparsity_pattern; + BlockSparseMatrix system_matrix; + + const unsigned int n_refinement_steps; + + double time_step; + double epsilon; + double vis; + double vfs_out; + double v_out; + + BlockVector solution; + BlockVector old_solution; + BlockVector 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 +class RightHandSide : public Function +{ + public: + RightHandSide () : Function(1) {}; + + virtual double value (const Point &p, + const unsigned int component = 0) const; +}; + + + +template +class PressureBoundaryValues : public Function +{ + public: + PressureBoundaryValues () : Function(1) {}; + + virtual double value (const Point &p, + const unsigned int component = 0) const; +}; + + +template +class SaturationBoundaryValues : public Function +{ + public: + SaturationBoundaryValues () : Function(dim+2) {}; + + virtual void vector_value (const Point &p, + Vector &value) const; +}; + + +template +class InitialValues : public Function +{ + public: + InitialValues () : Function(dim+2) {}; + + virtual double value (const Point &p, + const unsigned int component = 0) const; + + virtual void vector_value (const Point &p, + Vector &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 +double RightHandSide::value (const Point &/*p*/, + const unsigned int /*component*/) const +{ + return 0; +} + + + +template +double PressureBoundaryValues::value (const Point &p, + const unsigned int /*component*/) const +{ + return 1-p[0]; +} + + +template +void +SaturationBoundaryValues::vector_value (const Point &p, + Vector &values) const +{ + Assert (values.size() == dim+2, + ExcDimensionMismatch (values.size(), dim+2)); + + for (unsigned int d=0; d +double InitialValues::value (const Point &p, + const unsigned int component) const +{ + if(component +void +InitialValues::vector_value (const Point &p, + Vector &values) const +{ + Assert (values.size() == dim+2, + ExcDimensionMismatch (values.size(), dim+2)); + + for (unsigned int d=0; d +class KInverse : public TensorFunction<2,dim> +{ + public: + virtual void value_list (const std::vector > &points, + std::vector > &values, + const double epsilon) const; +}; + + +template +class Coefficient : public Function +{ + public: + Coefficient () : Function(1) {}; + + virtual double value (const Point &p, + const unsigned int component = 0 ) const; +}; + + +template +void +KInverse::value_list (const std::vector > &points, + std::vector > &values, + const double epsilon) const +{ + Assert (points.size() == values.size(), + ExcDimensionMismatch (points.size(), values.size())); + + for (unsigned int p=0; p +double Coefficient::value (const Point &/*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 +Tensor<1,dim> +extract_u (const FEValuesBase &fe_values, + const unsigned int i, + const unsigned int q) +{ + Tensor<1,dim> tmp; + + for (unsigned int d=0; d +double +extract_div_u (const FEValuesBase &fe_values, + const unsigned int i, + const unsigned int q) +{ + double divergence = 0; + for (unsigned int d=0; d +double extract_p (const FEValuesBase &fe_values, + const unsigned int i, + const unsigned int q) +{ + return fe_values.shape_value_component (i,q,dim); +} + +template +double extract_s (const FEValuesBase &fe_values, + const unsigned int i, + const unsigned int q) +{ + return fe_values.shape_value_component (i,q,dim+1); +} + +template +Tensor<1,dim> +extract_grad_s(const FEValuesBase &fe_values, + const unsigned int i, + const unsigned int q) +{ + Tensor<1,dim> tmp; + for (unsigned int d=0; d +TwoPhaseFlowProblem::TwoPhaseFlowProblem (const unsigned int degree) + : + degree (degree), + fe (FE_RaviartThomas(degree), 1, + FE_DGQ(degree), 1, + FE_DGQ(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 +void TwoPhaseFlowProblem::make_grid_and_dofs () +{ + GridGenerator::hyper_cube (triangulation, 0, 1); + + for (unsigned int f=0; f::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 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 +void TwoPhaseFlowProblem::assemble_system () +{ + QGauss quadrature_formula(degree+2); + QGauss face_quadrature_formula(degree+2); + + FEValues fe_values (fe, quadrature_formula, + update_values | update_gradients | + update_q_points | update_JxW_values); + FEFaceValues 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 local_matrix (dofs_per_cell, dofs_per_cell); + Vector local_rhs (dofs_per_cell); + + std::vector 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 right_hand_side; + const PressureBoundaryValues pressure_boundary_values; + const KInverse k_inverse; + const Coefficient coefficient; + + + std::vector rhs_values (n_q_points); + std::vector boundary_values (n_face_q_points); + std::vector > k_inverse_values (n_q_points); + std::vector coefficient_values(n_q_points); + + std::vector > old_solution_values(n_q_points, Vector(dim+2)); + std::vector > > old_solution_grads(n_q_points, + std::vector > (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::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 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 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::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 + 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 +class InverseMatrix : public Subscriptor +{ + public: + InverseMatrix (const Matrix &m); + + void vmult (Vector &dst, + const Vector &src) const; + + private: + const SmartPointer matrix; + + mutable GrowingVectorMemory<> vector_memory; +}; + + +template +InverseMatrix::InverseMatrix (const Matrix &m) + : + matrix (&m) +{} + + + +template +void InverseMatrix::vmult (Vector &dst, + const Vector &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 &A, + const InverseMatrix > &Minv); + + void vmult (Vector &dst, + const Vector &src) const; + + private: + const SmartPointer > system_matrix; + const SmartPointer > > m_inverse; + + mutable Vector tmp1, tmp2; +}; + + +SchurComplement::SchurComplement (const BlockSparseMatrix &A, + const InverseMatrix > &Minv) + : + system_matrix (&A), + m_inverse (&Minv), + tmp1 (A.block(0,0).m()), + tmp2 (A.block(0,0).m()) +{} + + +void SchurComplement::vmult (Vector &dst, + const Vector &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 &A); + + void vmult (Vector &dst, + const Vector &src) const; + + private: + const SmartPointer > system_matrix; + + mutable Vector tmp1, tmp2; +}; + + +ApproximateSchurComplement::ApproximateSchurComplement (const BlockSparseMatrix &A) + : + system_matrix (&A), + tmp1 (A.block(0,0).m()), + tmp2 (A.block(0,0).m()) +{} + + +void ApproximateSchurComplement::vmult (Vector &dst, + const Vector &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 +void TwoPhaseFlowProblem::solve () +{ + const InverseMatrix > + m_inverse (system_matrix.block(0,0)); + Vector tmp (solution.block(0).size()); + Vector schur_rhs (solution.block(1).size()); + Vector 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 + 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 quadrature_formula(degree+2); + QGauss face_quadrature_formula(degree+2); + FEValues fe_values (fe, quadrature_formula, + update_values | update_gradients | + update_q_points | update_JxW_values); + FEFaceValues fe_face_values (fe, face_quadrature_formula, + update_values | update_normal_vectors | + update_q_points | update_JxW_values); + FEFaceValues 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 local_rhs (dofs_per_cell); + std::vector > old_solution_values(n_q_points, Vector(dim+2)); + std::vector > old_solution_values_face(n_face_q_points, Vector(dim+2)); + std::vector > old_solution_values_face_neighbor(n_face_q_points, Vector(dim+2)); + std::vector > present_solution_values(n_q_points, Vector(dim+2)); + std::vector > present_solution_values_face(n_face_q_points, Vector(dim+2)); + + std::vector neighbor_saturation (n_face_q_points); + std::vector local_dof_indices (dofs_per_cell); + + typename DoFHandler::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 present_u; + for (unsigned int d=0; d 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::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::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;qat_boundary(face_no)) + { + if (cell->face(face_no)->boundary_indicator() ==2 ) + {for (unsigned int q=0;q present_u_face; + for (unsigned int d=0; d= 0); + + if (is_outflow_q_point == true) + { + for (unsigned int i=0; iget_dof_indices (local_dof_indices); + for (unsigned int i=0; i 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 +Vector +TwoPhaseFlowProblem::evaluate_solution (const Point &point) const +{ + static const MappingQ1 mapping; + // first find the cell in which this point + // is, initialize a quadrature rule with + // it, and then a FEValues object + const typename DoFHandler::active_cell_iterator + cell = GridTools::find_active_cell_around_point (dof_handler, point); + + const Point unit_point + = mapping.transform_real_to_unit_cell(cell, point); + Assert (GeometryInfo::is_inside_unit_cell (unit_point), + ExcInternalError()); + + const Quadrature quadrature (unit_point); + FEValues 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 > u_value(1, Vector(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 +void TwoPhaseFlowProblem::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 point(x,y); + + Vector solution_at_point(dim+2); + + solution_at_point = evaluate_solution (point); + + sampled_solution << point << " "; + for (unsigned int c=0; c +void TwoPhaseFlowProblem::output_results +(const unsigned int timestep_number) const +{ + std::vector 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 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 +void TwoPhaseFlowProblem::run () +{ + std::cout<<"Solving problem in " < production_rate; + std::list production_time; + + Vector tmp (old_solution.size()); + VectorTools::project (dof_handler, constraints, QGauss(degree+2),InitialValues(),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<::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; +}