opm-simulators/doc/handbook/tutorial1.tex
Andreas Lauser bf32eb0010 Implement the element centered finite volume spatial discretization
This makes eWoms multi-discretization capable. Along the way, this
fixes some bugs and does a medium sized reorganization of the source tree.

This is a squashed patch of the following commits:

--------
1st commit message:

add initial version of the element centered finite volume discretization

currently, it is a misnomer as it is just a copy of the vertex
centered discretization plus some renames...

--------
2nd commit message:

rename [VE]cfvModel -> [VE]cfvDiscretization

--------
3rd commit message:

ecfv: prelimary changes required to make it compile

but not work yet...

--------
4th commit message:

Rename *FvElementGeometry to *Stencil

"Stencil" seems to be the standard expression for this concept...

(also, it is not specific to finite volume methods and is shorter.)

--------
5th commit message:

refactor the stencil class for the element centered finite volume discretization

--------
6th commit message:

ECFV: some work on the stencil class

--------
7th commit message:

ECFV: make the boundary handling code compile

--------
8th commit message:

rename elemContext() to elementContext()

--------
9th commit message:

ECFV: make the VTK output modules compile

--------
10th commit message:

stencil: introduce the concept of primary DOFs

also save an vector of all element pointers in the stencil.

--------
11th commit message:

ECFV: try to fix assembly; add missing timeIdx arguments to the num*() methods

--------
12th commit message:

ECFV: fix stupid mistake in the assembler

--------
13th commit message:

ECFV: remove a few implicit DOF == vertex assumptions

the black-oil example now runs without valgrind complaints until it encounters
a negative oil mole fraction.

--------
14th commit message:

VCFV: make everything compile again

all vertex centered FV examples should now work again...

--------
15th commit message:

rename [ev]cfvmodel.hh to [ev]cfvdiscretization.hh

the classes have already been renamed.

--------
16th commit message:

ECFV: make it work to the point where it can write out the initial solution.

--------
17th commit message:

ECFV: make it work

the local residual/jacobian needed some work in distinguishing primary
and secondary DOFs and there was an minor issue with the serialization
code.

for some reason, it seems still not correct. (-> convergence is too slow.)

--------
18th commit message:

VCFV: make it compile for the black oil model again

--------
19th commit message:

VCFV: make it compile with the remaining models again

--------
20th commit message:

flash model: make it work with ECFV

although this breaks its compatibility with VCFV. (-> next commit)

--------
21st commit message:

adapt the VCFV to make it compatible with the flash model again

--------
22nd commit message:

make all models compile with VCFV again

--------
23rd commit message:

VCFV: more cleanups of the stencil

VcfvStencil now does not have any public attributes anymore. TODO: do
not export attributes in the SubControlVolume and SubControlVolumeFace
classes.

--------
24th commit message:

VCFV: actually update the element pointer

--------
25th commit message:

change the blackoil model back to ECFV

--------
26th commit message:

immiscible model: make it compatible with the ECFV discretization

--------
27th commit message:

PVS model: make it work with ECFV

--------
28th commit message:

NCP model: make it work with ECFV

--------
29th commit message:

rename Vcfv*VelocityModule to *VelocityModule

--------
30th commit message:

richards model: make it work with ECFV

--------
31st commit message:

unify the ECFV and the VCFV VTK output modules

and other cleanups

--------
32nd commit message:

unify the common code of the VCFV and the ECFV disctretizations

--------
33rd commit message:

unify the element contexts between element and vertex centered finite volumes

--------
34th commit message:

unify the local jacobian class of the finite volume discretizations

--------
35th commit message:

replace [VE]vcf(LocalResidual|ElementContext|BoundaryContext|ConstraintsContext) by generic code

--------
36th commit message:

replace the [EV]cfvLocalResidual by generic code

--------
37th commit message:

unify the MultiPhaseProblem and Problem classes, introduce NullBorderListCreator

--------
38th commit message:

remove the discretization specific boundary context

--------
39th commit message:

unify the [EV]cfvDiscretization classes

--------
40th commit message:

Unify [EV]cfvMultiPhaseFluxVariables

--------
41st commit message:

Unify the [EC]cfvNewton* classes

--------
42nd commit message:

Unify [EV]cfvVolumeVariables

--------
43rd commit message:

unify [EV]cfvAssembler

--------
44th commit message:

unified flux variables: fix stupid mistake when calculating pressure gradients

--------
45th commit message:

unify what's to unify for the [EV]CFV properties

--------
46th commit message:

make the method to calculate gradients and values at flux approximation points changeable

Currently, this is used by the vertex centered finite volume method to
be able to use P1-finite element gradients instead of two-point
ones...

--------
47th commit message:

make the restart code work correctly, use the correct DofMapper for VCFV

--------
48th commit message:

actually use the gradient calculator in a model

the immiscible model in this case

--------
49th commit message:

move some files around to where they belong, use the new gradient calculation code in all models

TODO: proper handling of boundary gradients

--------
50th commit message:

fix the stokes model

currently it only works with the vertex centered finite volume
discretization, but the plan is to soon move it to a staggered grid
scheme anyway...

--------
51st commit message:

move all models back to using the vertex centered finite volume discretization by default

--------
52nd commit message:

models: some variable renames and documentation fixes

- scv -> dof
- vert -> dof
- vertex -> dof
- replace 'VCFV'
- fix some typos

--------
53rd commit message:

don't expect UG anymore

since it is quite non-free and hard to get. we now use ALUGrid instead!

--------
54th commit message:

temporarily disable jacobian recycling

--------
55th commit message:

fix writing/reading restart files using the generic code

--------
56th commit message:

fix bug where fluxes were only counted once in the stencil

this only affected the vertex centered finite volumes discretization...

--------
57th commit message:

boundary gradients: use the center of the sub-control volume adjacent to a boundary segment

--------
58th commit message:

make it compile on GCC

--------
59th commit message:

get rid of most hacks

for this, partial reassemble and jacobian recycling was brought
back. For the this and the remaining stuff the main trick is the
introduction of the GridCommHandleFactory concept which constructs
communication handles suited for the respective spatial
discretization...

--------
60th commit message:

fix a few annoying bugs

first, default the convergence criterion for the linear solver did not
honor the initial residual which lead to linear solver breakdowns,
then some debugging code was left in the discrete fracture model and
then there was a bug in the TP gradient approximation class...

this has the consequence that we need a new reference solution for the
discrete fracture problem...

--------
61st commit message:

iterative linear solver: remove the code for the non-default convergence criteria

--------
62nd commit message:

provide the FE cache instead of the local FE

this fixes a segfault in the stokes model caused by the fact that the
local FE was not initialized at this point.

--------
63rd commit message:

(Navier-)Stokes: fix bug due to the transition to unit normals

now, all tests pass for this branch. The only things which need to be
fixed are some annoying performance regressions compared to master and
some bug in the splices feature of the property system...

--------
64th commit message:

some fix for the local residual of the immiscible model

--------
65th commit message:

Navier-Stokes: implement SCV center gradients

There seems to be a bug in the previous implementation (the jacobian
inverse transposed is evaluated using the local, not the global
geometry), so the reference solution for the stokes2c test problem has
also been updated...

--------
66th commit message:

remove the ALUGrid specialization of the LensGridCreator and the YaspGrid one for the fingerproblem

using different grid seems to sometimes cause a different vertex
order, which in turn causes the respective test to fail if the
reference solution was computed using the other grid...

--------
67th commit message:

VCFV: use the correct BorderListCreator

this makes MPI parallel computations work again. apart from
performance regressions, this branch does not exhibit any known
regressions compared to master anymore...

--------
68th commit message:

make verything compile with the element centered finite volume discretization

except the Navier-Stokes and the two-phase DFM models, of course...

--------
69th commit message:

minor fixes

- make the navier-stokes model slighly more generic by using the
  proper (in,ex)teriorIndex() methods on sub-control volumes
- make the signature of the calculateValue() template method of the
  common two-point gradient approximator match the one of the vertex
  centered finite volume one

--------
70th commit message:

fix fallout from the Big Rebase

--------
71st commit message:

ECFV: some bugs in the boundary

--------
72nd commit message:

make computeFlux() compute area-specific quantities

--------
73rd commit message:

fix more bugs in the element centered FV discretization

now eWoms should match Dumux pretty closely...

--------
74th commit message:

coalesce the common code of the multi phase porous medium models into "MultiPhaseBaseModel"

--------
75th commit message:

update reference solutions

these were changed because of the screw-up with the area of boundary
segments...

--------
76th commit message:

rename "ImplicitBase" to "FvBase"

because in eWoms, everything is implicit and these are currently the
base classes for all finite volume discretizations.

--------
77th commit message:

make the spatial discretization selectable using a splice

This requires an opm-core with a the patches from
https://github.com/OPM/opm-core/pull/446 merged...

--------
78th commit message:

rename the properties used for splices to *Splice

--------
79th commit message:

move the files in 'tests/models' to 'tests'

since 'tests' was empty except for the 'models' subdirectory...

--------
80th commit message:

improve and fix the tutorial

--------
81st commit message:

remove the -fno-strict-aliasing flag from the provided option files

seems like recent versions of Dune have been adapted...

--------
82nd commit message:

also compile all CO2 injection simulations using the element centered finite volume discretization

--------
83rd commit message:

PVS model: make it work properly with the element-centered finite volume discretiation

because DOF != number of vertices
2013-12-12 12:52:44 +01:00

549 lines
26 KiB
TeX

\section{Upfront considerations}
\label{tutorial1}
The process of setting up a problem using \eWoms can be roughly
divided into three parts:
\begin{enumerate}
\item A suitable model has to be chosen.
\item The geometry of the problem has to be defined and the suitable
grid has to be created.
\item Boundary and initial conditions and material properties have to
be specified.
\end{enumerate}
The physical set-up of the flow problem being solved in this tutorial
is illustrated in Figure \ref{tutorial1:problemfigure}: It
consists of a rectangular domain with no-flow boundary conditions on
the top and on the bottom of the domain, which is initially fully
saturated by a light oil. Water infiltrates from the left side and
replaces the oil. The effect of gravity is neglected.
\begin{figure}[ht]
\psfrag{x}{x}
\psfrag{y}{y}
\psfrag{no flow}{no flow}
\psfrag{water}{\textbf{water}}
\psfrag{oil}{\textcolor{white}{\textbf{oil}}}
\psfrag{p_w = 2 x 10^5 [Pa]}{$p_w = 2 \times 10^5$ [Pa]}
\psfrag{p_w_initial = 2 x 10^5 [Pa]}{\textcolor{white}{\textbf{$\mathbf{p_{w_{initial}} = 2 \times 10^5}$ [Pa]}}}
\psfrag{S_n = 0}{$S_n = 0$}
\psfrag{S_n_initial = 0}{\textcolor{white}{$\mathbf{S_{n_{initial}} = 1}$}}
\psfrag{q_w = 0 [kg/m^2s]}{$q_w = 0$ $\left[\frac{\textnormal{kg}}{\textnormal{m}^2 \textnormal{s}}\right]$}
\psfrag{q_n = -3 x 10^-4 [kg/m^2s]}{$q_n = 3 \times 10^{-2}$ $\left[\frac{\textnormal{kg}}{\textnormal{m}^2 \textnormal{s}}\right]$}
\centering
\includegraphics[width=0.9\linewidth,keepaspectratio]{EPS/tutorial-problemconfiguration}
\caption{Geometry of the tutorial problem with initial and boundary conditions.}\label{tutorial1:problemfigure}
\end{figure}
The solved equations are the mass balances of water and oil:
\begin{align}
\label{massbalancewater}
\frac {\partial (\phi \, S_{w}\, \varrho_{w})}{\partial t}
-
\nabla \cdot \left( \varrho_{w} \, \frac{k_{rw}}{\mu_{w}} \, \mathbf{K}\;\nabla p_w \right)
-
q_w
& =
0 \\
\label{massbalanceoil}
\frac {\partial (\phi \, S_{o}\, \varrho_{o})}{\partial t}
-
\nabla \cdot \left( \varrho_{o} \, \frac{k_{ro}}{\mu_{o}} \, \mathbf{K}\;\nabla p_o \right)
-
q_o
& =
0
\end{align}
\subsection{The main source file}
Listing \ref{tutorial1:mainfile} shows the main file
\texttt{tutorial/tutorial1.cc} for the tutorial problem using
a fully-implicit model that assumes immiscibility. This file has to be
compiled and executed in order to solve the problem described above.
\begin{lst}[File tutorial/tutorial1.cc]\label{tutorial1:mainfile} \mbox{}
\lstinputlisting[style=eWomsCode, numbersep=5pt, firstline=27, firstnumber=27]{../../tutorial/tutorial1.cc}
\end{lst}
In \eWoms, this file is usually quite short, as most of the work for
setting up the simulation is done by the generic startup routine
\texttt{Ewoms::start()}. The tasks left for the main source file are:
\begin{itemize}
\item Inclusion of the necessary header files from line
\ref{tutorial1:include-begin} to line
\ref{tutorial1:include-end}.
\item Specifying the type tag of the problem which is going to be
simulated at line \ref{tutorial1:set-type-tag}.
\item Starting the simulation using default \eWoms startup routine
\texttt{Ewoms::start()} on line \ref{tutorial1:call-start}.
\end{itemize}
The default \eWoms startup routine parses the command line arguments,
optionally reads a file which may specify further parameters, creates
the grid which represents the spatial domain, and then starts the
actual simulation. Parameters for the simulation, can be either
specified using command line arguments of the form
(\texttt{--parameter-name=value}), or in the file specified by the
\texttt{--parameter-file="file\_name"} argument. The list of all
parameters used by a simulation can be obtained by passing
\texttt{--help} to the executable on the command line.
\subsection{The problem file}
\begin{lst}[File tutorial/tutorial1problem.hh] \label{tutorial1:problemfile}\mbox{}
\lstinputlisting[style=eWomsCode, numbersep=5pt, firstline=28, firstnumber=28]{../../tutorial/tutorial1problem.hh}
\end{lst}
For using \eWoms, the central file is the \textit{problem file} as
shown in listing~\ref{tutorial1:problemfile}. This file is responsible
for specifying the physical setup of the problem which is to be
simulated. In this context, all problems first need to set up the
relevant properties for the \eWoms property system:
\begin{itemize}
\item First, a new type tag is created for the problem on line
\ref{tutorial1:create-type-tag}. In this case, the new type
tag inherits all properties from the \texttt{ImmiscibleTwoPhaseModel}
type tag, which means that for this problem the immiscible model
for two fluid phases is chosen as discretization scheme.
\item Next, the spatial discretization which ought to be used is
selected on line~\ref{tutorial1:include-discretization}. In this
case, the vertex-centered finite volume (VCFV) method is chosen.
\item On line \ref{tutorial1:set-problem}, the problem class is
specified for the new type tag, while the kind of grid that is to be
used is defined in line \ref{tutorial1:set-grid} -- in this case
it is \texttt{Dune::YaspGrid}.
\item Since \Dune does not provide a uniform mechanism to create or
load grids, \eWoms features the concept of grid creators. In this
case, the generic \texttt{CubeGridCreator} is used. This grid
creator constructs a structured rectangular grid with a specified
size and resolution. For this grid creator, the physical domain of
the grid is specified via the run-time parameters
\texttt{DomainSizeX} and \texttt{DomainSizeY} and its resolution by
\texttt{CellsX} and \texttt{CellsY}. These parameters can be
specified via the command-line or in a parameter file, but the
problem file must define default values for them. Defining these
defaults is done on lines \ref{tutorial1:grid-default-params-begin}
to \ref{tutorial1:default-params-end}.
\end{itemize}
Next, an appropriate fluid system -- which specifies the thermodynamic
relations of the fluid phases -- has to be chosen. By default, the
immiscible model for two fluid phases uses
\texttt{\justifyNoHyphen{}Ewoms::Fluid\-Systems::TwoPImmiscible}. This
fluid system simplifies things considerably by assuming immiscibility
of the phases, but it requires to explicitly specify the fluids of the
wetting and non-wetting phases. In this case, liquid water based on
simple relations is selected as the wetting phase on line
\ref{tutorial1:wettingPhase} and a light liquid oil is chosen as the
non-wetting phase on line \ref{tutorial1:nonwettingPhase}. The next
property, which is set on line \ref{tutorial1:gravity}, tells the
model not to use gravity. This is then followed by the specification
default values for parameters specifying the temporal and spatial
simulation domain from line \ref{tutorial1:default-params-begin} to
line \ref{tutorial1:default-params-end}. The values of those can be
overwritten at run-time if desired.
Following the property definitions, is the problem class which
specifies the variable parameters of the physical set-up to be
simulated. Such parameters are, for example, boundary and initial
conditions, source terms and spatially dependent quantities like
temperature, porosity, intrinsic permeability and the parameters
required by the capillary pressure/relative permeability law within
the domain. Since there are quite a few methods necessary to fully
specify a problem, and it is often possible to specify meaningful
defaults, it is strongly recomended to derive the problem from the
class specified by the \texttt{BaseProblem} property as done on line
\ref{tutorial1:def-problem}\footnote{Technically, deriving the problem
from the \texttt{BaseProblem} is not required if it implements all
necessary methods. Do not complain if things explode in this case,
though.}.
For the isothermal multi-phase porous media model which assumes
immiscibility, the problem class provided by the user must implement
at least the following methods:
\begin{itemize}
\item \texttt{initial()} for specifying the initial condition within
the domain,
\item \texttt{source()} which defines the source term for each conservation quantity,
\item \texttt{boundary()} for the conditions at the spatial domain's
boundary,
\item \texttt{temperature()} which provides the temperature within the
spatial domain. Note, that if energy is conserved, this method is
not necessary. Instead, energy fluxes must be specified by the
\texttt{boundary()} and \texttt{source()} methods and an initial
temperature needs to be defined by the \texttt{initial()} method in
this case.
\item \texttt{intrinsicPermeability()} for returning the intrinsic
permeability matrix $\mathbf{K}$ in $[m^2]$ at a given location.
\item \texttt{porosity()} which specifies the ratio of pore space to
total volume of the medium at a given location.
\item \texttt{materialLawParams()} which defines the parameters for
the capillary pressure and relative permeability relations at a
given location.
\end{itemize}
All of these methods take a single template argument,
\texttt{Context}, and the three function arguments \texttt{context},
\texttt{spaceIdx} and \texttt{timeIdx}. Together, these form the
so-called \textit{execution context}. The execution context can be
thought of as a collection of all available information for a given
method. Thus, execution contexts a way to abstract away the
differences of discretization schemes. The following methods are
provided by all \texttt{context} objects:
\begin{itemize}
\item \texttt{pos(spaceIdx, timeIdx)}: This method returns the
relevant position of the execution context within the spatial
domain.
\item \texttt{globalSpaceIndex(spaceIdx, timeIdx)}: Returns a global
index for the spatial degree of freedom which is associated to the
execution context. This index can be used to store spatially
dependent information in an array.
\item \texttt{element()}: Returns the \Dune grid element to which the
execution context applies.
\item \texttt{gridView()}: Returns the \Dune grid view of which the
element is part of.
\item \texttt{problem()}: Returns the \eWoms object which describes
the physical problem to be solved.
\item \texttt{model()}: Returns the \eWoms model which is used to simulate
the physical problem.
\end{itemize}
The execution contexts for the \texttt{source} and \texttt{boundary}
methods always provide the following methods in addition:
\begin{itemize}
\item \texttt{volVars(spaceIdx, timeIdx)}: The secondary variables
used by the model relevant for the execution context. These are
useful to specify solution dependent source and boundary terms.
\item \texttt{primaryVars(spaceIdx, timeIdx)}: The vector of primary
variables to which the model solves for.
\end{itemize}
Finally, the execution context for the \texttt{boundary} method
provides the methods \texttt{boundarySegmentArea(spaceIdx, timeIdx)},
and \texttt{normal(spaceIdx, timeIdx)}, which return the area of the
boundary segment and outer unit normal of the boundary segment.
When specifying the boundary conditions in the problem class, a small
safeguard value \texttt{eps\_} should usually be added when comparing
spatial coordinates. For the problem considered here, the left
boundary might not be detected if checking whether the first
coordinate of the global position is equal to zero, but it is detected
reliably by testing whether it is smaller than an
$\epsilon$-value. Also, the \texttt{bboxMax()} (``\textbf{max}imum
coordinate of the grid's \textbf{b}ounding \textbf{box}'') method is
used here to determine the extend of the physical domain. It returns a
vector with the maximum values on each axis of the grid. This method
and the analogous \texttt{bboxMin()} method are provided by the
problem's base class.
In addition to these methods, there might be some additional
model-specific methods.
\subsection{Defining fluid properties}
\label{tutorial1:description-fluid-class}
The \eWoms distribution includes representations of some common
substances which can be used out of the box. The properties of pure
substances (such as nitrogen, water, or the pseudo-component air) are
provided by header files located in the folder
\texttt{ewoms/material/components}.
When two or more substances are considered, interactions between these
become relevant and the properties of these mixtures of substances --
e.g. composition, density or enthalpy -- are required. In \eWoms, such
interactions are defined by {\em fluid systems}, which are located in
\texttt{ewoms/material/fluidsystems}. A more thorough overview of the
\eWoms fluid framework can be found in
chapter~\ref{sec:fluidframework}.
\subsection{Exercises}
\label{tutorial1:exercises}
The following exercises will give you the opportunity to learn how you
can change soil parameters, boundary conditions, run-time parameters
and fluid properties in \eWoms.
\subsubsection{Exercise 1}
\renewcommand{\labelenumi}{\alph{enumi})}
For Exercise 1 you have to modify the tutorial files only sightly, or
even not at all.
\begin{enumerate}
\item \textbf{Running the Program} \\
To get an impression what the results should look like, you can
first run the original version of the tutorial problem by entering
the directory into which you installed \eWoms, followed typing
\texttt{make tutorial1 \&\& ./bin/tutorial1}. Note, that the
time-step size is automatically adapted during the simulation. For
visualizing the results using the program \texttt{paraview}, please
refer to section \ref{quick-start-guide}.
\item \textbf{Changing the Model Domain} \\
Change the size of the model domain so that you get a rectangle with
edge lengths of $\text{x} = \unit[400]{m}$ and $\text{y} =
\unit[500]{m}$ and with discretization lengths of $\Delta \text{x} =
\unit[20]{m}$ and $\Delta \text{y} = \unit[20]{m}$. For this, you can
either change the properties defined between lines
\ref{tutorial1:default-params-begin} and
\ref{tutorial1:default-params-end}, or you may pass them as
command line parameters to the simulation when you run it.
If you chose to change the problem file, re-compile the simulation
by typing \texttt{make tutorial1} in the toplevel directory of your
\eWoms build directory. Note, that you do not have to recompile the
program if you use command line parameters. In both cases, don't
forget run the simulation before re-inspecting the results using
\texttt{paraview}.
\item \textbf{Boundary Conditions} \\
Change the boundary conditions in the problem file
\texttt{tutorial1problem.hh} so that water enters from the
bottom and oil is extracted from the top boundary. The right and the
left boundary should be closed for both, water and oil.
Note that in \eWoms, flux boundary conditions are specified as
fluxes of the conserved quantities into the direction of the outer
normal per area. This means a mass flux into the domain is negative,
a flux out of the domain is always positive. Re-compile the
simulation by typing \texttt{make tutorial1} and re-run it
as explained above.
\item \textbf{Changing the Spatial Discretization} \\
Change the spatial discretization used for this problem from the
vertex-centered finite element method (VCFV) to the element-centered
one (ECFV). For this, you need to change the file included on
line~\ref{tutorial1:include-discretization} and the definition of
the splice for the spatial discretization on
line~\ref{tutorial1:set-spatial-discretization}.
\item \textbf{Changing the Shape of the Discrete Elements} \\
In order to complete this exercise you need a \Dune grid manager
that is capable of handling simplex grids. By default, \Dune does
not include such a grid manager in its core modules, but the freely
available \texttt{ALUGrid}~\cite{ALUGRID-HP} can be used. If you do
not want to go through the trouble of installing this grid manager,
please skip this exercise.
Change the grid creator used by the problem to
\texttt{SimplexGridCreator<TypeTag>} and the type of the grid to
\texttt{Dune::ALUSimplexGrid<2, 2>}. The grid creator is specified
on line \ref{tutorial1:set-gridcreator}, whil the type of the
\Dune grid manager is set on line
\ref{tutorial1:set-grid}. You also need to change the include
statement of the grid creator from \texttt{cubegridcreator.hh} to
\texttt{simplexgridcreator.hh} on line
\ref{tutorial1:include-grid-creator} and the one for the grid
manager from \texttt{dune/grid/yaspgrid.hh} to
\texttt{dune/grid/alugrid.hh} on line \ref{tutorial1:include-grid-manager}.
The resulting grid can be examined by re-compiling and starting the
simulation, loading the result into \texttt{paraview}, and selecting
\texttt{Surface with Edges} instead of the default visualization
mode \texttt{Surface}.
\item \textbf{Changing Fluids} \\
In this exercise you will change the fluids used: Use DNAPL instead
of LNAPL and Brine instead of Water. This can be done by the
following steps:
\begin{enumerate}
\item Brine: Brine is thermodynamically very similar to pure water but
also includes a fixed amount of salt in the liquid phase. Hence,
the class \texttt{Ewoms::Brine} calls back to a class which
represents pure water. (This can be, for example
\texttt{Ewoms::H2O}, or \texttt{Ewoms::SimpleH2O}.) The class which
represents pure water is passed to \texttt{Ewoms::Brine} as the
second template argument after the data type for scalar values,
i.e. the full definition of the brine component is
\texttt{Ewoms::Brine<Scalar, Ewoms::H2O<Scalar> >}. The file which
defines the brine component is located in the folder
\texttt{ewoms/material/components/}. Include this file and select
use the brine component as the wetting phase.
\item DNAPL: Now let us use a component that represents an oil that
is heavier than water (in environmental engineering this class of
substances is often called \textbf{d}ense
\textbf{n}on-\textbf{a}queous \textbf{p}hase \textbf{l}iquid
(DNAPL)) which is also located in the folder
\texttt{ewoms/material/components/}. Include the correct file and
use the DNAPL component as the non-wetting phase.
\end{enumerate}
If you want to take a closer look on how the fluid classes are defined
and which substances are available in the \eWoms distribution, feel
free to browse through the files in the directory
\texttt{ewoms/material/components} and read
chapter~\ref{sec:fluidframework}.
\item \textbf{Use a Full-Fledged Fluid System} \\
In \eWoms, the canonical way to describe fluid mixtures is via
\textit{fluid systems}\footnote{For a thorough introduction into
fluid systems and the concepts related to it, see chapter
\ref{sec:fluidframework}}. In order to include a fluid system,
you first have to comment out lines
\ref{tutorial1:2p-system-start}
to \ref{tutorial1:2p-system-end} in the problem file.\\
Now include the file
\texttt{ewoms/material/fluidsystems/h2oairfluidsystem.hh}, and
instruct the model to use this fluid system by setting
the \texttt{FluidSystem} property via:\\
\begin{lstlisting}[style=eWomsCode]
SET_TYPE_PROP(Tutorial1Problem,
FluidSystem,
Ewoms::FluidSystems::H2OAir<typename GET_PROP_TYPE(TypeTag, Scalar)>)
\end{lstlisting}
However, this is a rather complicated fluid system which considers all
fluid phases as a mixture of the components and also uses tabulated
components that need to be filled with values (i.e. some components
use tables to speed up expensive computations). The initialization of
the fluid system is normally done in the constructor of the problem by
calling \texttt{FluidSystem::init();}.
As water flow replacing a gas is much faster, simulate the problem
only until $2.000$ seconds is reached and start with a time step of
$1$ second.
Now, revert the changes made in this part of the exercise, as we will
continue to use immiscible phases from here on and hence do not need a
complex fluid system.
\item \textbf{Changing Constitutive Relations} \\
Instead of using a regularized Brooks-Corey law for the relative
permeability and for the capillary pressure saturation relationship,
use an linear material law which does not consider residual
saturation. The new material law should exhibit an entry pressure of
$p_e = 0.0\;\text{Pa}$ and maximal capillary pressure of
$p_{c_{max}} = 2000.0\;\text{Pa}$. To do that, you have to change
the material law property on line
\ref{tutorial1:eff2abs} as follows:
\begin{itemize}
\item First, switch to the linear material law. This can be achieved
by replacing \texttt{RegularizedBrooksCorey} in the private section
of the property definition by \texttt{LinearMaterial} and rename
\texttt{RawMaterialLaw} to \texttt{TwoPMaterialLaw}. Also remove the
line which contains the \texttt{EffToAbsLaw} typedef.
\item Then, adapt the necessary parameters of the linear law and the
respective \texttt{set}-functions can be found in the file
\texttt{ewoms/material/fluidmatrixinteractions/2p/linearmaterialparams.hh}.\\
Call the \texttt{set}-functions from the constructor of the problem.
\item Finally, re-compile and re-run the simulation. In order to
successfully compile the simulation, you also have to include the
header file which contains the linear material law. This header is
called \texttt{linearmaterial.hh} and is located in the same
directory as the header for the regularized Brooks-Corey law,
\texttt{regularizedbrookscorey.hh}.
\end{itemize}
\item \textbf{Heterogeneities} \\
Set up a model domain with the soil properties given in Figure
\ref{tutorial1:exercise1_d}. Adjust the boundary conditions
so that water is flowing from the left to the right of the domain again.
\begin{figure}[ht]
\psfrag{K1 =}{$\mathbf{K} = 10^{-8}\;\text{m}^2$}
\psfrag{phi1 =}{$\phi = 0.15$}
\psfrag{K2 =}{\textcolor{white}{$\mathbf{K} = 10^{-9}\;\text{m}^2$}}
\psfrag{phi2 =}{\textcolor{white}{$\phi = 0.3$}}
\psfrag{600 m}{$600 \;\text{m}$}
\psfrag{300 m}{$300 \;\text{m}$}
\centering
\includegraphics[width=0.5\linewidth,keepaspectratio]{EPS/exercise1_c.eps}
\caption{Exercise 1f: Set-up of a model domain with a heterogeneity. $\Delta x = 20 \;\text{m}$ $\Delta y = 20\;\text{m}$.}\label{tutorial1:exercise1_d}
\end{figure}
You can use the fluids of exercise 1b).
\textbf{Hint:} The relevant position in the domain can be obtained using
\texttt{const auto \&pos=context.pos(spaceIdx, timeIdx);}
When does the front cross the material border? In paraview, the
animation view (\textit{View} $\rightarrow$ \textit{Animation View})
is a convenient way to get a rough feeling of the time-step sizes.
\end{enumerate}
\subsubsection{Exercise 2}
For this exercise, first create a new problem file analogous to the
file \texttt{tutorial1problem.hh} (e.g. with the name
\texttt{ex2\_tutorial1problem.hh}. The new problem
file needs to be included in the file \texttt{tutorial1.cc}.
The new file should contain definitions of new classes with names that
relate to the file name, such as
\texttt{Ex2TutorialProblemCoupled}. Make sure that you also adjust the
guardian macros in lines \ref{tutorial1:guardian1} and
\ref{tutorial1:guardian1} in the header files (e.g. change
\mbox{\texttt{EWOMS\_TUTORIAL1PROBLEM\_HH}} to
\mbox{\texttt{EWOMS\_EX2\_TUTORIAL1PROBLEM\_HH}}). Include the
new problem file in \texttt{tutorial1.cc}. Besides adjusting
the guardian macros, the new problem file should define and use a new
type tag for the problem as well as a new problem class
e.g. \mbox{\texttt{Ex2TutorialProblemCoupled}}. The type tag for the
problem should be adjusted, too. For this, modify line
\ref{tutorial1:set-type-tag} in the problem file and the adapt
the \texttt{main} function in the file \texttt{tutorial1.cc}.
After this, change the domain parameters so that they match the domain
described by figure \ref{tutorial1:ex2_Domain}. Adapt the
problem class so that the boundary conditions are consistent with
figure \ref{tutorial1:ex2_BC}. Initially, the domain is fully
saturated with water and the pressure is $p_w = 5\cdot
10^5\;\text{Pa}$. Oil infiltrates from the left side. Create a grid
with $20$ cells in $x$-direction and $10$ cells in $y$-direction. The
simulation time should end at $10^6\;\text{s}$ with an initial time
step size of $100\;\text{s}$.
\begin{figure}[ht]
\psfrag{K1}{K $= 10^{-7}\;\text{m}^2$}
\psfrag{phi1}{$\phi = 0.2$}
\psfrag{Lin}{\textsc{Brooks}-\textsc{Corey} Law}
\psfrag{Lin2}{$\lambda = 1.8$, $p_e = 1000\;\text{Pa}$}
\psfrag{K2}{K $= 10^{-9}\;\text{m}^2$}
\psfrag{phi2}{$\phi = 0.15$}
\psfrag{BC1}{\textsc{Brooks}-\textsc{Corey} Law}
\psfrag{BC2}{$\lambda = 2$, $p_e = 1500\;\text{Pa}$}
\psfrag{H1y}{$50\;\text{m}$}
\psfrag{H2y}{$15\;\text{m}$}
\psfrag{H3y}{$20\;\text{m}$}
\psfrag{L1x}{$100\;\text{m}$}
\psfrag{L2x}{$50\;\text{m}$}
\psfrag{L3x}{$25\;\text{m}$}
\centering
\includegraphics[width=0.8\linewidth,keepaspectratio]{EPS/Ex2_Domain.eps}
\caption{Set-up of the model domain and the soil parameters}\label{tutorial1:ex2_Domain}
\end{figure}
\begin{figure}[ht]
\psfrag{pw}{$p_w = 5 \times 10^5\;\text{Pa}$}
\psfrag{S}{$S_n = 1.0$}
\psfrag{qw}{$q_w = 2 \times 10^{-4} \;\text{kg}/\text{m}^2\text{s}$}
\psfrag{qo}{$q_n = 0.0 \;\text{kg}/\text{m}^2\text{s}$}
\psfrag{no flow}{no flow}
\centering
\includegraphics[width=0.8\linewidth,keepaspectratio]{EPS/Ex2_Boundary.eps}
\caption{Boundary Conditions}\label{tutorial1:ex2_BC}
\end{figure}
\begin{itemize}
\item Increase the simulation time to e.g. $4\times 10^7
\;\text{s}$. Investigate the saturation of the wetting phas: Is the
value range reasonable?
\item What happens if you use a grid with more cells in each direction?
\end{itemize}
\subsubsection{Exercise 3: Create a New Component}
Create a new file for the benzene component called \texttt{benzene.hh}
and implement a new component. (\textbf{Hint:} The easiest way to do
this is to copy and modify one of the existing components located in
the directory \texttt{ewoms/material/components}.) Use benzene as a
new non-wetting fluid and run the problem of Exercise 2 with water and
benzene. Benzene has a density of $889.51 \, \text{kg} / \text{m}^3$
and a viscosity of $0.00112 \, \text{Pa} \; \text{s}$.
\clearpage \newpage
%%% Local Variables:
%%% mode: latex
%%% TeX-master: "ewoms-handbook"
%%% End: