tutorial-coupled.tex 27.8 KB
Newer Older
1
\section[Fully-Implicit Model]{Solving a Problem Using a Fully-Coupled Model}\label{tutorial-coupled}
Bernd Flemisch's avatar
Bernd Flemisch committed
2

3
The process of setting up a problem using \Dumux can be roughly divided into four parts:
Bernd Flemisch's avatar
Bernd Flemisch committed
4
\begin{enumerate}
5
 \item A suitable model has to be chosen.
Bernd Flemisch's avatar
Bernd Flemisch committed
6
 \item The geometry of the problem and correspondingly a grid have to be defined.
7
 \item Material properties and constitutive relationships have to be selected.
8
 \item Boundary conditions and initial conditions have to be specified.
Bernd Flemisch's avatar
Bernd Flemisch committed
9
10
\end{enumerate}

11
The problem being solved in this tutorial is illustrated in Figure \ref{tutorial-coupled:problemfigure}. 
12
A rectangular domain with no-flow boundaries on the top and on the bottom, which is initially saturated with oil, is considered. 
13
Water infiltrates from the left side into the domain and replaces the oil. Gravity effects are neglected here.
Bernd Flemisch's avatar
Bernd Flemisch committed
14

15
\begin{figure}[ht]
Bernd Flemisch's avatar
Bernd Flemisch committed
16
17
18
19
20
21
22
23
24
25
\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]$}
26
\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]$}
Bernd Flemisch's avatar
Bernd Flemisch committed
27
28
29
30
31
\centering
\includegraphics[width=0.9\linewidth,keepaspectratio]{EPS/tutorial-problemconfiguration}
\caption{Geometry of the tutorial problem with initial and boundary conditions.}\label{tutorial-coupled:problemfigure}
\end{figure}

32
The solved equations are the mass balances of water and oil:
Bernd Flemisch's avatar
Bernd Flemisch committed
33
34
35
\begin{align}
  \label{massbalancewater}
  \frac {\partial (\phi \, S_{w}\, \varrho_{w})}{\partial t}
36
  -
37
  \nabla \cdot \left( \varrho_{w} \, \frac{k_{rw}}{\mu_{w}} \, \mathbf{K}\;\nabla p_w \right)
Bernd Flemisch's avatar
Bernd Flemisch committed
38
  -
39
  q_w
Bernd Flemisch's avatar
Bernd Flemisch committed
40
41
42
43
  & =
  0 \\
  \label{massbalanceoil}
  \frac {\partial (\phi \, S_{o}\, \varrho_{o})}{\partial t}
44
  -
45
  \nabla \cdot \left( \varrho_{o} \, \frac{k_{ro}}{\mu_{o}} \, \mathbf{K}\;\nabla p_o \right)
Bernd Flemisch's avatar
Bernd Flemisch committed
46
  -
47
  q_o 
Bernd Flemisch's avatar
Bernd Flemisch committed
48
49
50
51
  & =
  0
\end{align}

52
\subsection{The Main File}
Bernd Flemisch's avatar
Bernd Flemisch committed
53

54
Listing \ref{tutorial-coupled:mainfile} shows the main application file
55
\texttt{tutorial/tutorial\_coupled.cc} for the coupled two-phase
56
model. This file has to be compiled and executed in order to solve the problem described
Bernd Flemisch's avatar
Bernd Flemisch committed
57
58
59
above.

\begin{lst}[File tutorial/tutorial\_coupled.cc]\label{tutorial-coupled:mainfile} \mbox{}
Christoph Grueninger's avatar
Christoph Grueninger committed
60
  \lstinputlisting[style=DumuxCode, numbersep=5pt, firstline=29]{../../tutorial/tutorial_coupled.cc}
Bernd Flemisch's avatar
Bernd Flemisch committed
61
62
63
\end{lst}

From line \ref{tutorial-coupled:include-begin} to line
64
\ref{tutorial-coupled:include-end} the required headers are included.
Bernd Flemisch's avatar
Bernd Flemisch committed
65
66

At line \ref{tutorial-coupled:set-type-tag} the type tag of the
67
68
problem, which is going to be simulated, is specified. All other data
types can be retrieved via the \Dumux property system and only depend
69
on this single type tag. For a more thorough introduction to the
70
71
\Dumux property system, see chapter~\ref{sec:propertysystem}.

72
After this, the default startup routine \texttt{Dumux::start()} is
73
74
called on line \ref{tutorial-coupled:call-start}. This function deals
with parsing the command line arguments, reading the parameter file,
75
76
77
78
79
setting up the infrastructure necessary for \Dune, loading the grid, and
starting the simulation. 
Required parameters for the start of the simulation, 
such as the initial time-step size, the simulation time or details of the grid,
can be either specified by command line arguments of the form
80
81
(\texttt{-ParameterName ParameterValue}), in the file specified by the
\texttt{-parameterFile} argument, or if the latter is not specified,
82
83
in the file \mbox{\texttt{tutorial\_coupled.input}}. 
If a parameter is
84
85
86
87
88
89
specified on the command line as well as in the parameter file, the
values provided in the command line have
precedence. Listing~\ref{tutorial-coupled:parameter-file} shows the
default parameter file for the tutorial problem.

\begin{lst}[File tutorial/tutorial\_coupled.input]\label{tutorial-coupled:parameter-file} \mbox{}
Christoph Grueninger's avatar
Christoph Grueninger committed
90
\lstinputlisting[style=DumuxCode]{../../tutorial/tutorial_coupled.input}
91
92
93
94
95
\end{lst}

To provide an error message, the usage message which is displayed to
the user if the simulation is called incorrectly, is printed via the
custom function which is defined on
96
97
98
99
line~\ref{tutorial-coupled:usage-function}. 
In this function the usage message is customized to the problem at hand. 
This means that at least the necessary parameters are listed here.  
For more information about the input file please refer to section \ref{sec:inputFiles}.
100

101
102

\subsection{The Problem Class}
Bernd Flemisch's avatar
Bernd Flemisch committed
103
104

When solving a problem using \Dumux, the most important file is the
105
106
so-called \textit{problem file} as shown in
listing~\ref{tutorial-coupled:problemfile}.
Bernd Flemisch's avatar
Bernd Flemisch committed
107
108

\begin{lst}[File tutorial/tutorialproblem\_coupled.hh]\label{tutorial-coupled:problemfile} \mbox{}
Christoph Grueninger's avatar
Christoph Grueninger committed
109
\lstinputlisting[style=DumuxCode, numbersep=5pt, firstline=28]{../../tutorial/tutorialproblem_coupled.hh}
Bernd Flemisch's avatar
Bernd Flemisch committed
110
111
\end{lst}

112
First, a new type tag is created for the problem in line
Bernd Flemisch's avatar
Bernd Flemisch committed
113
\ref{tutorial-coupled:create-type-tag}.  In this case, the new type
114
115
116
117
118
119
120
121
122
123
124
125
tag inherits all properties from the \texttt{BoxTwoP} type tag, which
means that for this problem the two-phase box model is chosen as
discretization scheme. Further, it inherits from the spatial
parameters type tag, which is defined in the problem-dependent spatial
parameters file (line
\ref{tutorial-coupled:define-spatialparameters-typetag}).  On line
\ref{tutorial-coupled:set-problem}, a problem class is attached to the
new type tag, while the grid which is going to be used is defined in
line \ref{tutorial-coupled:set-grid} -- in this case that is
\texttt{Dune::YaspGrid}.  Since there's no uniform mechanism to
allocate grids in \Dune, \Dumux features the concept of grid creators.
In this case the generic \texttt{CubeGridCreator} which creates a
126
127
structured hexahedron grid of a specified size and resolution. For
this grid creator the  physical domain of the grid is specified via the
128
129
130
run-time parameters \texttt{Grid.upperRightX},
\texttt{Grid.upperRightY}, \texttt{Grid.numberOfCellsX} and
\texttt{Grid.numberOfCellsY}. These parameters can be specified via
131
the command-line or in a parameter file.
132
133
134
135

Next, the appropriate fluid system, which specifies the thermodynamic
relations of the fluid phases, has to be chosen. By default, the
two-phase model uses the \texttt{TwoPImmiscibleFluidSystem}, which
136
assumes immiscibility of the phases, but requires the components
137
138
used for the wetting and non-wetting phases to be explicitly set. In
this case, liquid water which uses the relations from
139
IAPWS'97~\cite{IAPWS1997} is chosen as the wetting phase on line
140
141
142
143
144
145
146
147
148
149
150
151
\ref{tutorial-coupled:wettingPhase} and liquid oil is chosen as the
non-wetting phase on line \ref{tutorial-coupled:nonwettingPhase}. The
last property, which is set in line \ref{tutorial-coupled:gravity},
tells the model not to use gravity.

Parameters which are specific to a physical set-up to be simulated,
such as boundary and initial conditions, source terms or temperature
within the domain, and which are required to solve the differential
equations of the models are specified via a \textit{problem} class. If
the two-phase box model is used, this class should be derived from
\texttt{TwoPBoxProblem} as done in line
\ref{tutorial-coupled:def-problem}.
Bernd Flemisch's avatar
Bernd Flemisch committed
152
153
154

The problem class always has at least five methods:
\begin{itemize}
155
156
\item A method \texttt{boundaryTypes()} specifying the type of
  boundary conditions at each vertex.
Bernd Flemisch's avatar
Bernd Flemisch committed
157
\item A method \texttt{dirichlet()} specifying the actual values for
158
  the \textsc{Dirichlet} conditions at each \textsc{Dirichlet} vertex.
Bernd Flemisch's avatar
Bernd Flemisch committed
159
\item A method \texttt{neumann()} specifying the actual values for
160
161
  the \textsc{Neumann} conditions, which are usually evaluated at the 
  integration points of the \textsc{Neumann} boundary faces.
162
\item A method for source or sink terms called \texttt{source()}, usually evaluated at
163
  the center of a control volume.
Bernd Flemisch's avatar
Bernd Flemisch committed
164
\item A method called \texttt{initial()} for specifying the initial
165
  conditions at each vertex.
Bernd Flemisch's avatar
Bernd Flemisch committed
166
167
\end{itemize}

168
169
170
For the definition of the the boundary condition types and of the
values of the \textsc{Dirichlet} boundaries, two parameters are
available:
171
172
\begin{description}
 \item [values:]  A vector which stores the result of the method. What
173
  the values in this vector mean is dependent on the method: For
174
175
  \texttt{dirichlet()} it contains the actual values of the primary
  variables, for \texttt{boundaryTypes()} it contains the boundary 
176
177
178
  condition types. It has as many entries as the model has primary variables / equations.
  For the typical case, in which all equations have the same boundary
  condition type at a certain position, there are two methods that set the appropriate conditions
179
180
181
182
183
  for all primary variables / equations: \texttt{setAllDirichlet()} and \texttt{setAllNeumann()}.
\item [vertex:] The boundary condition and the Dirichlet values are
  specified for a vertex, which represents a control volume in the box
  discretization. This avoids the specification of two different
  boundary condition types for one equation at different control
184
  volumes.  Be aware that the second parameter is a Dune grid entity
185
  with the codimension \texttt{dim}.
186
187
\end{description}

188
189
To ensure that no boundaries are undefined, a small safeguard value
\texttt{eps\_} is usually added when comparing spatial
190
191
coordinates. The left boundary is hence not detected by checking, if the
first coordinate of the global position is equal to zero, but by testing whether it is
192
193
smaller than a very small value \texttt{eps\_}.

194
Methods which make statements about boundary segments of the grid
195
(such as \texttt{neumann()}) are called with six arguments:
Bernd Flemisch's avatar
Bernd Flemisch committed
196
\begin{description}
197
198
\item[values:] A vector \texttt{neumann()}, in which the mass fluxes per area unit
  over the boundary segment are specified.
Bernd Flemisch's avatar
Bernd Flemisch committed
199
200
201
202
\item[element:] The element of the grid where the boundary segment
  is located.
\item[fvElemGeometry:] The finite-volume geometry induced on the
  finite element by the box scheme.
203
\item[isIt:] The \texttt{Intersection} of the boundary segment as given by the grid.
Bernd Flemisch's avatar
Bernd Flemisch committed
204
\item[scvIdx:] The index of the sub-control volume in
205
  \texttt{fvElementGeometry} which is assigned to the boundary segment.
Bernd Flemisch's avatar
Bernd Flemisch committed
206
207
208
209
\item[boundaryFaceIdx:] The index of the boundary face in
  \texttt{fvElementGeometry} which represents the boundary segment.  
\end{description}

Melanie Darcis's avatar
Melanie Darcis committed
210
Similarly, the \texttt{initial()} and \texttt{source()} methods
211
specify properties of control volumes and thus only get
Bernd Flemisch's avatar
Bernd Flemisch committed
212
\texttt{values}, \texttt{element}, \texttt{fvElemGeom} and
213
\texttt{scvIdx} as arguments.
Bernd Flemisch's avatar
Bernd Flemisch committed
214
215

In addition to these five methods, there might be some model-specific
216
217
218
219
220
221
methods. If the isothermal two-phase model is used, this includes for
example a \texttt{temperature()} method which returns the temperature
in \textsc{Kelvin} of the fluids and the rock matrix in the
domain. This temperature is then used by the model to calculate fluid
properties which possibly depend on it, e.g. density. The
\texttt{bboxMax()} (``\textbf{max}imum coordinated of the grid's
222
223
224
\textbf{b}ounding \textbf{b}ox'') method is used here to
determine the extend of the physical domain. It returns a vector with the
maximum values of each global coordinate of the grid. This method
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
and the analogous \texttt{bboxMin()} method are provided by the base
class \texttt{Dumux::BoxProblem<TypeTag>}.

\subsection{Defining Fluid Properties}\label{tutorial-coupled:description-fluid-class}

The \Dumux distribution includes some common substances which can be
used out of the box. The properties of the pure substances (such as
the components nitrogen, water, or the pseudo-component air) are
provided by header files located in the folder
\verb+dumux/material/components+.

Most often, when two or more components are considered, fluid
interactions such as solubility effects come into play and properties
of mixtures such as density or enthalpy are of interest. These
interactions are defined by {\em fluid systems}, which are located in
\verb+dumux/material/fluidsystems+. A more thorough overview of the
\Dumux fluid framework can be found in chapter~\ref{sec:fluidframework}.
Bernd Flemisch's avatar
Bernd Flemisch committed
242

243
244
% In this example, a class for the definition of a two-phase system is used. This allows for the choice 
% of the two components oil and water and for access of the parameters that are relevant for the two-phase model.
245

246
\subsection{Defining Spatially Dependent Parameters}\label{tutorial-coupled:description-spatialParameters}
Bernd Flemisch's avatar
Bernd Flemisch committed
247

248
249
250
251
252
In \Dumux, many properties of the porous medium can depend on the
spatial location. Such properties are the \textit{intrinsic
  permeability}, the parameters of the \textit{capillary pressure} and
the \textit{relative permeability}, the \textit{porosity}, the
\textit{heat capacity} as well as the \textit{heat conductivity}. Such
253
parameters are defined using a so-called \textit{spatial parameters}
254
class.
Bernd Flemisch's avatar
Bernd Flemisch committed
255

256
If the box discretization is to used, the spatial parameters class
257
should be derived from the base class
258
\texttt{Dumux::BoxSpatialParams<TypeTag>}. Listing
259
\ref{tutorial-coupled:spatialparametersfile} shows the file \\
260
\verb+tutorialspatialparams_coupled.hh+:
Bernd Flemisch's avatar
Bernd Flemisch committed
261

262
263
\begin{lst}[File tutorial/tutorialspatialparams\_coupled.hh]\label{tutorial-coupled:spatialparametersfile} \mbox{}
\lstinputlisting[style=DumuxCode, numbersep=5pt, firstline=28]{../../tutorial/tutorialspatialparams_coupled.hh}
Bernd Flemisch's avatar
Bernd Flemisch committed
264
265
\end{lst}

266
267
First, the spatial parameters type tag is created on line
\ref{tutorial-coupled:define-spatialparameters-typetag}. The type tag
268
269
for the problem is then derived from it. The \Dumux properties defined on
the type tag for the spatial parameters are, for example, the spatial
270
271
parameters class itself (line
\ref{tutorial-coupled:set-spatialparameters}) or the capillary
272
pressure/relative permeability relations\footnote{Taken together, the
273
274
275
276
277
278
279
280
  capillary pressure and the relative permeability relations are
  called \textit{material law}.} which ought to be used by the
simulation (line
\ref{tutorial-coupled:rawlaw} \label{tutorial-coupled:materialLaw}).
\Dumux provides several material laws in the folder
\verb+dumux/material/fluidmatrixinteractions+.  The selected one --
here it is a relation according to a regularized version of
\textsc{Brooks} \& \textsc{Corey} -- is included in line
281
282
283
284
285
\ref{tutorial-coupled:rawLawInclude}. 
After the selection, an adapter class is specified in line \ref{tutorial-coupled:eff2abs} to
translate between effective and absolute saturations. Like this,
residual saturations can be specified in a generic way.  As only the employed
material law knows the names of the parameters which it
286
requires, it provides a parameter class
287
288
\texttt{RegularizedBrooksCoreyParams} which has the type
\texttt{Params} and which is defined in line
289
290
\ref{tutorial-coupled:matLawObjectType}. In this case, the spatial
parameters only require a single set of parameters which means that it
291
only requires a single material parameter object as can be seen in
292
293
294
295
line~\ref{tutorial-coupled:matParamsObject}.

In line \ref{tutorial-coupled:permeability}, a method returning the
intrinsic permeability is specified. As can be seen, the method has
296
to be called with three arguments:
297
298
299
\begin{description}
\item[\texttt{element}:] Just like for the problem itself, this
  parameter describes the considered element by means of a \Dune
300
  entity. Elements provide information about their geometry and
301
  position and can be mapped to a global index.
302
303
304
305
306
\item[\texttt{fvElemGeom}:] It holds information about the finite-volume
  geometry of the element induced by the box method.
\item[\texttt{scvIdx}:] This is the index of the sub-control volume of the
  element which is considered. It is equivalent to the local index
  of the vertex which corresponds to the considered control volume in
307
308
309
  the element.
\end{description}

310
311
The intrinsic permeability is usually a tensor. Thus the method returns
a $\texttt{dim} \times \texttt{dim}$-matrix, where \texttt{dim} is the
312
dimension of the grid.
313

314
The method \texttt{porosity()} defined in line
Bernd Flemisch's avatar
Bernd Flemisch committed
315
\ref{tutorial-coupled:porosity} is called with the same arguments as
316
317
318
\texttt{intrinsicPermeability()} and returns a scalar value for
porosity dependent on the position in the domain.

319
320
321
Next, the method \texttt{materialLawParams()}, defined in line
\ref{tutorial-coupled:matLawParams}, returns the
\verb+materialLawParams+ object that is applied at the specified
322
323
324
325
326
327
328
329
330
331
332
333
position. Although in this case only one object is returned, in
general, the problem may be heterogeneous, which necessitates
returning different objects at different positions in space.  While
the selection of the type of this object was already explained (line
\ref{tutorial-coupled:rawLawInclude}), some specific parameter values
of the used material law, such as the \textsc{Brooks} \&
\textsc{Corey} parameters, are still needed. This is done in the
constructor at line \ref{tutorial-coupled:setLawParams}.  Depending on
the type of the \texttt{materialLaw} object, the \texttt{set}-methods
might be different than those given in this example. The name of the
access / set functions as well as the rest of the implementation of
the material description can be found in
334
\verb+dumux/material/fluidmatrixinteractions/2p+.
Bernd Flemisch's avatar
Bernd Flemisch committed
335
336
337
338

\subsection{Exercises}
\label{tutorial-coupled:exercises}
The following exercises will give you the opportunity to learn how you
339
340
can change soil parameters, boundary conditions, run-time parameters
and fluid properties in \Dumux.
Bernd Flemisch's avatar
Bernd Flemisch committed
341
342

\subsubsection{Exercise 1}
343
344
\renewcommand{\labelenumi}{\alph{enumi})} For Exercise 1 you have
to make only some small changes in the tutorial files.  
345

Bernd Flemisch's avatar
Bernd Flemisch committed
346
347
\begin{enumerate}

348
\item \textbf{Running the Program} \\
Christoph Grueninger's avatar
Christoph Grueninger committed
349
  To get an impression what the results should look like you can first run the original version of the coupled tutorial model by typing  \texttt{./tutorial\_coupled}. Note, that the time-step size is automatically adapted during the simulation. For the visualization of the results using paraview please refer to section \ref{quick-start-guide}.
350
351

\item \textbf{Changing the Model Domain and the Boundary Conditions} \\
Bernd Flemisch's avatar
Bernd Flemisch committed
352
  Change the size of the model domain so that you get a rectangle with
353
354
355
  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}$.
Bernd Flemisch's avatar
Bernd Flemisch committed
356
357
358
359
360
361
362
  
  Change the boundary conditions in the file
  \texttt{tutorialproblem\_coupled.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 water and oil fluxes. 

  Compile the main file by typing \texttt{make tutorial\_coupled} and
363
  run the model as explained above.
Bernd Flemisch's avatar
Bernd Flemisch committed
364

365
  \item \textbf{Changing  the Shape of the Discrete Elements} \\
366
367
  Change the types of elements used for discretizing the domain. In line \ref{tutorial-coupled:set-gridcreator} of the problem file  the type of gridcreator is chosen. By choosing a different grid creator you can discretize the domain with different elements. Hint: You can find gridcreators in \texttt{dumux/common/}. The shape of the employed elements can be visualized in paraview by choosing \texttt{Surface with Edges}. 

Bernd Flemisch's avatar
Bernd Flemisch committed
368
\item \textbf{Changing Fluids} \\
369
Now you can change the fluids. Use DNAPL instead of Oil and Brine instead of Water. To do that, you have to select different components via the property system in the problem file:
370
\begin{enumerate}
371
 \item Brine: Brine is thermodynamically very similar to pure water but also considers a fixed amount of salt in the liquid phase. Hence, the class \texttt{Dumux::Brine} uses a pure water class, such as \texttt{Dumux::H2O}, as a second template argument after the data type \texttt{<Scalar>}.
372
 \item DNAPL: A standard set of chemical substances, such as Oil and Brine, is already included in the problem (via a list of \texttt{\#include ..} statements) and hence easily accessible by default. However, this is not the case for the class \texttt{Dumux::SimpleDNAPL}, which describes a simple \textbf{d}ense \textbf{n}on-\textbf{a}queous \textbf{p}hase \textbf{l}iquid and is located in the folder \texttt{dumux/material/components/}. Try to include the file and select the component as the non-wetting phase via the property system.
373
\end{enumerate}
374
If you want to take a closer look on how the fluid classes are defined and which substances are already available please browse through the files in the directory
375
\texttt{/dumux/material/components} and read chapter~\ref{sec:fluidframework}.
Bernd Flemisch's avatar
Bernd Flemisch committed
376

377
\item \textbf{Use a Full-Fledged Fluid System} \\
378
\Dumux usually describes fluid mixtures via \textit{fluid systems}, see also chapter \ref{sec:fluidframework}. In order to include a fluid system, you first have to comment out lines \ref{tutorial-coupled:2p-system-start} to \ref{tutorial-coupled:2p-system-end} in the problem file. If you use eclipse, this can easily be done by pressing \textit{Ctrl + Shift + 7} -- the same as to cancel the comment later on.\\
379
Now include the file \texttt{fluidsystems/h2oairfluidsystem.hh} in the material folder, and set a property \texttt{FluidSystem} with the appropriate type, i.e. \texttt{Dumux::H2OAirFluidSystem<TypeTag>}. However, this is a rather complicated fluid system which considers mixtures of components and also uses tabulated components that need to be initialized -- i.e. the tables need to be filled with values. The initialization of the fluid system is normally done in the constructor of the problem by calling \texttt{GET\_PROP\_TYPE(TypeTag, FluidSystem)::init();}. As water flow replacing a gas is much faster, test your simulation only until $2000$ seconds and start with a time step of $1$ second.\\
380
Please reverse 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.
Bernd Flemisch's avatar
Bernd Flemisch committed
381
382

\item \textbf{Changing Constitutive Relations} \\
383
  Use an unregularized linear law with an entry pressure of $p_e = 0.0\;\text{Pa}$ and maximal capillary pressure of e.g. $p_{c_{max}} = 2000.0\;\text{Pa}$ instead of using a
384
 regularized Brooks-Corey law for the
385
  relative permeability and for the capillary pressure saturation relationship. To do that you have
386
  to change the file \texttt{tutorialspatialparams\_coupled.hh}. 
Bernd Flemisch's avatar
Bernd Flemisch committed
387
 You can find the material laws in the folder 
Melanie Darcis's avatar
Melanie Darcis committed
388
  \verb+dumux/material/fluidmatrixinteractions+. The necessary parameters
389
390
391
of the linear law and the respective \texttt{set}-functions can be found
 in the file \\
 \verb+dumux/material/fluidmatrixinteractions/2p/linearmaterialparams.hh+.
Bernd Flemisch's avatar
Bernd Flemisch committed
392
393
394
395
 
\item \textbf{Heterogeneities}  \\
  Set up a model domain with the soil properties given in Figure
  \ref{tutorial-coupled:exercise1_d}. Adjust the boundary conditions
396
  so that water is again flowing from the left to the right of the
397
\begin{figure}[ht]
398
\psfrag{K1 =}{$\mathbf{K} = 10^{-8}\;\text{m}^2$}
Bernd Flemisch's avatar
Bernd Flemisch committed
399
\psfrag{phi1 =}{$\phi = 0.15$}
400
\psfrag{K2 =}{\textcolor{white}{$\mathbf{K} = 10^{-9}\;\text{m}^2$}}
Bernd Flemisch's avatar
Bernd Flemisch committed
401
\psfrag{phi2 =}{\textcolor{white}{$\phi = 0.3$}}
402
403
\psfrag{600 m}{$600 \;\text{m}$}
\psfrag{300 m}{$300 \;\text{m}$}
Bernd Flemisch's avatar
Bernd Flemisch committed
404
405
\centering
\includegraphics[width=0.5\linewidth,keepaspectratio]{EPS/exercise1_c.eps}
406
\caption{Exercise 1f: Set-up of a model domain with a heterogeneity. $\Delta x = 20 \;\text{m}$ $\Delta y = 20\;\text{m}$.}\label{tutorial-coupled:exercise1_d}
Bernd Flemisch's avatar
Bernd Flemisch committed
407
\end{figure}
408
domain. You can use the fluids of exercise 1c).\\
409
\textbf{Hint:} The current position of the control volume can be obtained using \texttt{element\allowbreak.geometry()\allowbreak.corner(scvIdx)}.\\
410
When does the front cross the material border? In paraview, the
411
animation view (\textit{View} $\rightarrow$ \textit{Animation
412
413
  View}) is a convenient way to get a rough feeling of the time-step
sizes.
Bernd Flemisch's avatar
Bernd Flemisch committed
414
415
416
\end{enumerate}

\subsubsection{Exercise 2}
417
For this exercise you should create a new problem file analogous to
418
419
the file \texttt{tutorialproblem\_coupled.hh} (e.g. with the name
\texttt{ex2\_tutorialproblem\_coupled.hh} and new spatial parameters
420
just like \texttt{tutorialspatialparams\_coupled.hh}. The new
421
problem file needs to
422
423
be included in the file \texttt{tutorial\_coupled.cc}.

424
425
426
427
428
The new files should contain the definition 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{tutorial-coupled:guardian1} and
\ref{tutorial-coupled:guardian1}
429
430
431
432
in the header files (e.g. change
\mbox{\texttt{DUMUX\_TUTORIALPROBLEM\_COUPLED\_HH}} to\\
\mbox{\texttt{DUMUX\_EX2\_TUTORIALPROBLEM\_COUPLED\_HH}}). 
Besides adjusting the guardian macros, the new problem file should define and
433
use a new type tag for the problem as well as a new problem class
434
e.g. \mbox{\texttt{Ex2TutorialProblemCoupled}}. Make sure to assign your
435
newly defined spatial parameter class to the
436
\texttt{SpatialParams} property for the new
437
438
type tag. 

439
440
441
442
443
444
445
446
447
After this, change the run-time parameters so that they match the
domain described by figure \ref{tutorial-coupled:ex2_Domain}. Adapt
the problem class so that the boundary conditions are consistent with
figure \ref{tutorial-coupled:ex2_BC}. Initially, the domain is fully
saturated with water and the pressure is $p_w = 5 \times
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 be set to $10^6\;\text{s}$ with an
initial time step size of $100\;\text{s}$.
Bernd Flemisch's avatar
Bernd Flemisch committed
448
449
450
451
452
453

Now include your new problem file in the main file and replace the
\texttt{TutorialProblemCoupled} type tag by the one you've created and
compile the program.


454
\begin{figure}[ht]
455
\psfrag{K1}{K $= 10^{-7}\;\text{m}^2$}
Bernd Flemisch's avatar
Bernd Flemisch committed
456
\psfrag{phi1}{$\phi = 0.2$}
457
\psfrag{Lin}{\textsc{Brooks}-\textsc{Corey} Law}
458
459
\psfrag{Lin2}{$\lambda = 1.8$, $p_e = 1000\;\text{Pa}$}
\psfrag{K2}{K $= 10^{-9}\;\text{m}^2$}
Bernd Flemisch's avatar
Bernd Flemisch committed
460
\psfrag{phi2}{$\phi = 0.15$}
461
462
463
464
465
466
467
468
\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}$}
Bernd Flemisch's avatar
Bernd Flemisch committed
469
470
471
472
473
\centering
\includegraphics[width=0.8\linewidth,keepaspectratio]{EPS/Ex2_Domain.eps}
\caption{Set-up of the model domain and the soil parameters}\label{tutorial-coupled:ex2_Domain}
\end{figure}

474
\begin{figure}[ht]
475
\psfrag{pw}{$p_w = 5 \times 10^5\;\text{Pa}$}
Bernd Flemisch's avatar
Bernd Flemisch committed
476
\psfrag{S}{$S_n = 1.0$}
477
478
\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}$}
Bernd Flemisch's avatar
Bernd Flemisch committed
479
480
481
482
483
484
\psfrag{no flow}{no flow}
\centering
\includegraphics[width=0.8\linewidth,keepaspectratio]{EPS/Ex2_Boundary.eps}
\caption{Boundary Conditions}\label{tutorial-coupled:ex2_BC}
\end{figure}

485
\begin{itemize}
486
 \item Increase the simulation time to e.g. $4\times 10^7 \;\text{s}$. Investigate the saturation: Is the value range reasonable?
487
488
489
 \item What happens if you increase the resolution of the grid?
\end{itemize}

490
\subsubsection{Exercise 3: Parameter File Input}
491
492
493

As you have experienced, compilation takes quite some time. Therefore,
\Dumux provides a simple method to read in parameters at run-time
Christoph Grueninger's avatar
Christoph Grueninger committed
494
via \textit{parameter input files}.
495

496
497
498
499
500
501
502
In the code, parameters can be read via the macro
\texttt{GET\_RUNTIME\_PARAM(TypeTag, Scalar,
  MyWonderfulGroup.MyWonderfulParameter);}. At the end of the
simulation a list of parameters is printed if the command line option
\texttt{-printParams 1} is passed to the simulation. Add some (for
example \texttt{Newton.MaxSteps} and \texttt{EnableGravity}) to the
parameter file \texttt{tutorial\_coupled.input} and observe what
503
happens if they are modified. For more information about the input file please refer to section \ref{sec:inputFiles}.
504

505
\subsubsection{Exercise 4: Create a New Component}
Bernd Flemisch's avatar
Bernd Flemisch committed
506

507
508
Create a new file for the benzene component called \texttt{benzene.hh}
and implement a new component. (You may get a hint by looking at
509
existing components in the directory \verb+/dumux/material/components+). \\
Bernd Flemisch's avatar
Bernd Flemisch committed
510
Use benzene as a new fluid and run the model of Exercise 2 with water
511
512
and benzene. Benzene has a density of $889.51 \, \text{kg} /
\text{m}^3$ and a viscosity of $0.00112 \, \text{Pa} \; \text{s}$.
Bernd Flemisch's avatar
Bernd Flemisch committed
513

514
\clearpage \newpage
515
516
517
518
%%% Local Variables: 
%%% mode: latex
%%% TeX-master: "dumux-handbook"
%%% End: