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

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

Thomas Fetzer's avatar
Thomas Fetzer committed
12
13
The problem being solved in this tutorial is illustrated in Figure \ref{tutorial-coupled:problemfigure}.
A rectangular domain with no-flow boundaries on the top and on the bottom, which is initially saturated with oil, is considered.
14
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
15

16
\begin{figure}[ht]
Bernd Flemisch's avatar
Bernd Flemisch committed
17
\centering
18
19
\begin{tikzpicture}[>=latex]
  % basic sketch
Thomas Fetzer's avatar
Thomas Fetzer committed
20
  \fill [fill=dumuxBlue](0,0) rectangle ++(2,1.5);
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
  \fill [fill=dumuxYellow](2,0) rectangle ++(5,1.5);
  \draw (0,0) rectangle ++(7,1.5);
  \foreach \x in {0,0.25,...,6.75}
    \draw (\x,1.5) -- ++(0.25,0.25);
  \foreach \x in {0,0.25,...,6.75}
    \draw (\x,-0.25) -- ++(0.25,0.25);
  % labels
  \draw [->](-0.5,-0.5) -- ++(0,0.7) node [anchor=east]{$y$};
  \draw [->](-0.5,-0.5) -- ++(0.7,0) node [anchor=north]{$x$};
  \node at(3.5,1.75)[anchor=south]{no flow};
  \node at(3.5,-0.25)[anchor=north]{no flow};
  \draw[->,thick](-0.4,0.75) -- ++(0.9,0) node[anchor=north]{water};
  \draw[->,thick](6.5,0.75)node[anchor=north]{oil} -- ++(0.9,0);
  % equations
  \node [anchor=west] at (2,1.1){$p_{w_\text{initial}} = \unit[2 \cdot 10^5]{Pa}$};
  \node [anchor=west] at (2,0.4){$S_{n_\text{initial}} = 1$};
  \node [anchor=west] at (-3,1.1){$p_w = \unit[2 \cdot 10^5]{Pa}$};
  \node [anchor=west] at (-3,0.4){$S_n = 0$};
  \node [anchor=west] at (7.5,1.1){$q_w = \unitfrac[0]{kg}{m^2s}$};
  \node [anchor=west] at (7.5,0.4){$q_n = \unitfrac[3 \cdot 10^{-2}] {kg}{m^2s}$};
\end{tikzpicture}
\caption{Geometry of the tutorial problem with initial and boundary conditions.}
\label{tutorial-coupled:problemfigure}
Bernd Flemisch's avatar
Bernd Flemisch committed
44
45
\end{figure}

46
The solved equations are the mass balances of water and oil:
Bernd Flemisch's avatar
Bernd Flemisch committed
47
48
49
\begin{align}
  \label{massbalancewater}
  \frac {\partial (\phi \, S_{w}\, \varrho_{w})}{\partial t}
50
  -
51
  \nabla \cdot \left( \varrho_{w} \, \frac{k_{rw}}{\mu_{w}} \, \mathbf{K}\;\nabla p_w \right)
Bernd Flemisch's avatar
Bernd Flemisch committed
52
  -
53
  q_w
Bernd Flemisch's avatar
Bernd Flemisch committed
54
55
56
57
  & =
  0 \\
  \label{massbalanceoil}
  \frac {\partial (\phi \, S_{o}\, \varrho_{o})}{\partial t}
58
  -
59
  \nabla \cdot \left( \varrho_{o} \, \frac{k_{ro}}{\mu_{o}} \, \mathbf{K}\;\nabla p_o \right)
Bernd Flemisch's avatar
Bernd Flemisch committed
60
  -
Thomas Fetzer's avatar
Thomas Fetzer committed
61
  q_o
Bernd Flemisch's avatar
Bernd Flemisch committed
62
63
64
65
  & =
  0
\end{align}

66
\subsection{The Main File}
67
Listing \ref{tutorial-coupled:mainfile} shows the main application file
68
\texttt{tutorial/tutorial\_coupled.cc} for the coupled two-phase
69
model. This file has to be compiled and executed in order to solve the problem described
Bernd Flemisch's avatar
Bernd Flemisch committed
70
71
72
above.

\begin{lst}[File tutorial/tutorial\_coupled.cc]\label{tutorial-coupled:mainfile} \mbox{}
73
  \lstinputlisting[style=DumuxCode, numbersep=5pt, firstline=24, firstnumber=24]{../../tutorial/tutorial_coupled.cc}
Bernd Flemisch's avatar
Bernd Flemisch committed
74
75
76
\end{lst}

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

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

85
After this, the default startup routine \texttt{Dumux::start()} is
86
87
called on line \ref{tutorial-coupled:call-start}. This function deals
with parsing the command line arguments, reading the parameter file,
88
setting up the infrastructure necessary for \Dune, loading the grid, and
Thomas Fetzer's avatar
Thomas Fetzer committed
89
90
starting the simulation.
Required parameters for the start of the simulation,
91
92
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
93
(\texttt{-ParameterName ParameterValue}), in the file specified by the
94
\texttt{-ParameterFile} argument, or if the latter is not specified,
Thomas Fetzer's avatar
Thomas Fetzer committed
95
in the file \mbox{\texttt{tutorial\_coupled.input}}.
96
If a parameter is
97
98
99
100
101
102
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{}
Thomas Fetzer's avatar
Thomas Fetzer committed
103
\lstinputlisting[style=DumuxParameterFile]{../../tutorial/tutorial_coupled.input}
104
105
106
107
108
\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
Thomas Fetzer's avatar
Thomas Fetzer committed
109
110
111
line~\ref{tutorial-coupled:usage-function} in the main file.
In this function the usage message is customized to the problem at hand.
This means that at least the necessary parameters are listed here.
Bernd Flemisch's avatar
Bernd Flemisch committed
112

113
114
\subsection{The Problem Class}
\label{tutorial-coupled:problem}
Bernd Flemisch's avatar
Bernd Flemisch committed
115
When solving a problem using \Dumux, the most important file is the
116
117
so-called \textit{problem file} as shown in
listing~\ref{tutorial-coupled:problemfile}.
Bernd Flemisch's avatar
Bernd Flemisch committed
118
119

\begin{lst}[File tutorial/tutorialproblem\_coupled.hh]\label{tutorial-coupled:problemfile} \mbox{}
120
\lstinputlisting[style=DumuxCode, numbersep=5pt, firstline=24, firstnumber=24]{../../tutorial/tutorialproblem_coupled.hh}
Bernd Flemisch's avatar
Bernd Flemisch committed
121
122
\end{lst}

123
First, a new type tag is created for the problem in line
Bernd Flemisch's avatar
Bernd Flemisch committed
124
\ref{tutorial-coupled:create-type-tag}.  In this case, the new type
125
126
127
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
128
129
130
parameters type tag, which is defined in line
\ref{tutorial-coupled:define-spatialparameters-typetag} of the problem-dependent spatial
parameters file.  On line
131
132
133
134
135
136
\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
137
138
structured hexahedron grid of a specified size and resolution. For
this grid creator the  physical domain of the grid is specified via the
139
140
141
run-time parameters \texttt{Grid.upperRightX},
\texttt{Grid.upperRightY}, \texttt{Grid.numberOfCellsX} and
\texttt{Grid.numberOfCellsY}. These parameters can be specified via
142
the command-line or in a parameter file.
143
144
145
146

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
147
assumes immiscibility of the phases, but requires the components
148
149
used for the wetting and non-wetting phases to be explicitly set. In
this case, liquid water which uses the relations from
150
IAPWS'97~\cite{IAPWS1997} is chosen as the wetting phase on line
151
152
153
154
155
\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.

156
\label{tutorial-coupled:boundaryStart}Parameters which are specific to a physical set-up to be simulated,
157
158
such as boundary and initial conditions, source terms or temperature
within the domain, and which are required to solve the differential
Thomas Fetzer's avatar
Thomas Fetzer committed
159
equations of the models are specified via a \textit{problem} class. This class
160
should be derived from \texttt{ImplicitPorousMediaProblem} as done in line
161
\ref{tutorial-coupled:def-problem}.
Bernd Flemisch's avatar
Bernd Flemisch committed
162
163
164

The problem class always has at least five methods:
\begin{itemize}
165
166
\item A method \texttt{boundaryTypes()} specifying the type of
  boundary conditions at each vertex.
Bernd Flemisch's avatar
Bernd Flemisch committed
167
\item A method \texttt{dirichlet()} specifying the actual values for
168
  the \textsc{Dirichlet} conditions at each \textsc{Dirichlet} vertex.
Bernd Flemisch's avatar
Bernd Flemisch committed
169
\item A method \texttt{neumann()} specifying the actual values for
Thomas Fetzer's avatar
Thomas Fetzer committed
170
  the \textsc{Neumann} conditions, which are usually evaluated at the
171
  integration points of the \textsc{Neumann} boundary faces.
172
\item A method for source or sink terms called \texttt{source()}, usually evaluated at
173
  the center of a control volume.
Bernd Flemisch's avatar
Bernd Flemisch committed
174
\item A method called \texttt{initial()} for specifying the initial
175
  conditions at each vertex.
Bernd Flemisch's avatar
Bernd Flemisch committed
176
177
\end{itemize}

178
For the definition of the boundary condition types and of the
179
180
values of the \textsc{Dirichlet} boundaries, two parameters are
available:
181
\begin{description}
182
 \item [bcTypes/values:]  A vector which stores the result of the method. What
183
  the values in this vector mean is dependent on the method: For
184
  \texttt{dirichlet()}, \texttt{values} contains the actual values of the primary
Thomas Fetzer's avatar
Thomas Fetzer committed
185
  variables, for \texttt{boundaryTypes()}, \texttt{bcTypes} contains the boundary
186
187
188
  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
189
  for all primary variables / equations: \texttt{setAllDirichlet()} and \texttt{setAllNeumann()}.
190
191
192
193
194
\item [vertex:] The boundary condition and the Dirichlet values are
  specified for a vertex, which represents a sub-control volume in the box
  discretization. This inhibits the specification of two different
  boundary condition types for one equation at one sub-control
  volume.  Be aware that the second parameter is a Dune grid entity
195
  with the co-dimension \texttt{dim}.
196
197
\end{description}

198
199
To ensure that no boundaries are undefined, a small safeguard value
\texttt{eps\_} is usually added when comparing spatial
200
201
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
202
203
smaller than a very small value \texttt{eps\_}.

204
Methods for box models which make statements about boundary segments of the grid
205
(such as \texttt{neumann()}) are called with six arguments:
Bernd Flemisch's avatar
Bernd Flemisch committed
206
\begin{description}
207
208
\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
209
210
\item[element:] The element of the grid where the boundary segment
  is located.
211
\item[fvGeometry:] The finite-volume geometry induced on the
Bernd Flemisch's avatar
Bernd Flemisch committed
212
  finite element by the box scheme.
213
\item[intersection:] The \texttt{Intersection} of the boundary segment as given by the grid.
Bernd Flemisch's avatar
Bernd Flemisch committed
214
\item[scvIdx:] The index of the sub-control volume in
215
  \texttt{fvGeometry} which is assigned to the boundary segment.
Bernd Flemisch's avatar
Bernd Flemisch committed
216
\item[boundaryFaceIdx:] The index of the boundary face in
Thomas Fetzer's avatar
Thomas Fetzer committed
217
  \texttt{fvGeometry} which represents the boundary segment.
Bernd Flemisch's avatar
Bernd Flemisch committed
218
219
\end{description}

Melanie Darcis's avatar
Melanie Darcis committed
220
Similarly, the \texttt{initial()} and \texttt{source()} methods
221
specify properties of control volumes and thus only get
222
\texttt{values}, \texttt{element}, \texttt{fvGeometry} and
223
\texttt{scvIdx} as arguments.
Bernd Flemisch's avatar
Bernd Flemisch committed
224
225

In addition to these five methods, there might be some model-specific
226
227
228
229
230
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
231
232
\texttt{bBoxMax()} (``\textit{max}imum coordinated of the grid's
\textit{b}ounding \textit{box}'') method is used here to
233
234
determine the extend of the physical domain. It returns a vector with the
maximum values of each global coordinate of the grid. This method
235
and the analogous \texttt{bBoxMin()} method are provided by the base
236
237
class \texttt{Dumux::BoxProblem<TypeTag>}.

238
239
\subsection{Defining Fluid Properties}
\label{tutorial-coupled:description-fluid-class}
240
241
242
243
244
245
246
247
248
249
250
251
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
252

Thomas Fetzer's avatar
Thomas Fetzer committed
253
% In this example, a class for the definition of a two-phase system is used. This allows for the choice
254
% of the two components oil and water and for access of the parameters that are relevant for the two-phase model.
255

256
257
\subsection{Defining Spatially Dependent Parameters}
\label{tutorial-coupled:description-spatialParameters}
258
259
260
261
262
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
263
parameters are defined using a so-called \textit{spatial parameters}
264
class.
Bernd Flemisch's avatar
Bernd Flemisch committed
265

266
If the box discretization is used, the spatial parameters class
267
should be derived from the base class
268
\texttt{Dumux::BoxSpatialParams<TypeTag>}. Listing
269
\ref{tutorial-coupled:spatialparametersfile} shows the file \\
270
\verb+tutorialspatialparams_coupled.hh+:
Bernd Flemisch's avatar
Bernd Flemisch committed
271

272
\begin{lst}[File tutorial/tutorialspatialparams\_coupled.hh]\label{tutorial-coupled:spatialparametersfile} \mbox{}
273
\lstinputlisting[style=DumuxCode, numbersep=5pt, firstline=25, firstnumber=25]{../../tutorial/tutorialspatialparams_coupled.hh}
Bernd Flemisch's avatar
Bernd Flemisch committed
274
275
\end{lst}

276
277
First, the spatial parameters type tag is created on line
\ref{tutorial-coupled:define-spatialparameters-typetag}. The type tag
278
279
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
280
281
parameters class itself (line
\ref{tutorial-coupled:set-spatialparameters}) or the capillary
282
pressure/relative permeability relations\footnote{Taken together, the
283
284
285
286
287
288
289
290
  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
Thomas Fetzer's avatar
Thomas Fetzer committed
291
\ref{tutorial-coupled:rawLawInclude}.
292
293
294
295
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
296
requires, it provides a parameter class
297
298
\texttt{RegularizedBrooksCoreyParams} which has the type
\texttt{Params} and which is defined in line
299
300
\ref{tutorial-coupled:matLawObjectType}. In this case, the spatial
parameters only require a single set of parameters which means that it
301
only requires a single material parameter object as can be seen in
302
303
304
305
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
306
to be called with three arguments:
307
308
309
\begin{description}
\item[\texttt{element}:] Just like for the problem itself, this
  parameter describes the considered element by means of a \Dune
310
  entity. Elements provide information about their geometry and
311
  position and can be mapped to a global index.
312
\item[\texttt{fvGeometry}:] It holds information about the finite-volume
313
314
315
316
  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
317
318
319
  the element.
\end{description}

320
321
The intrinsic permeability is usually a tensor. Thus the method returns
a $\texttt{dim} \times \texttt{dim}$-matrix, where \texttt{dim} is the
322
dimension of the grid.
323

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

329
330
Next, the method \texttt{materialLawParams()}, defined in line
\ref{tutorial-coupled:matLawParams}, returns the
331
\verb+materialParams_+ object that is applied at the specified
332
333
334
335
336
337
338
339
340
341
342
343
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
344
\verb+dumux/material/fluidmatrixinteractions/2p+.
Bernd Flemisch's avatar
Bernd Flemisch committed
345
346
347
348

\subsection{Exercises}
\label{tutorial-coupled:exercises}
The following exercises will give you the opportunity to learn how you
349
can change soil parameters, boundary conditions, run-time parameters
350
and fluid properties in \Dumux. Possible solutions to these exercises are given in the tutorial folder in the
Thomas Fetzer's avatar
Thomas Fetzer committed
351
352
sub-folder \texttt{solutions\_coupled} as \texttt{.diff} files. In these files only
the lines that are different from the original file are listed.
353
They can be opened using the program \texttt{kompare}, simply type \texttt{kompare SOLUTIONFILE} into the terminal.
Bernd Flemisch's avatar
Bernd Flemisch committed
354
355

\subsubsection{Exercise 1}
356
\renewcommand{\labelenumi}{\alph{enumi})} For Exercise 1 you have
Thomas Fetzer's avatar
Thomas Fetzer committed
357
to make only some small changes in the tutorial files.
358

Bernd Flemisch's avatar
Bernd Flemisch committed
359
360
\begin{enumerate}

361
\item \textbf{Running the Program} \\
Thomas Fetzer's avatar
Thomas Fetzer committed
362
To get an impression what the results should look like you can compile and run the original version of
363
the coupled tutorial model by typing \texttt{make tutorial\_coupled} followed by \texttt{./tutorial\_coupled}.
Thomas Fetzer's avatar
Thomas Fetzer committed
364
Note, that the time-step size is automatically adapted during the simulation.
365
For the visualization of the results using ParaView, please refer to section \ref{quick-start-guide}.
366
367

\item \textbf{Changing the Model Domain and the Boundary Conditions} \\
Bernd Flemisch's avatar
Bernd Flemisch committed
368
  Change the size of the model domain so that you get a rectangle with
369
370
  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
Thomas Fetzer's avatar
Thomas Fetzer committed
371
  \text{y} = \unit[20]{m}$. For this you have to edit the parameter file (\texttt{tutorialproblem\_coupled.input})
372
373
  and run the program again.\\
  Note, that you do not have to recompile the program if you make changes to the parameter file.
374

Thomas Fetzer's avatar
Thomas Fetzer committed
375

Bernd Flemisch's avatar
Bernd Flemisch committed
376
377
378
  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
379
  left boundary should be closed for water and oil fluxes. \\
Thomas Fetzer's avatar
Thomas Fetzer committed
380
381
  The Neumannn Boundary conditions are multiplied by the normal (pointing outwards),
  so an influx is negative, an outflux always positive.
382
  Such information can easily be found in the documentation of the functions (also look into base classes).
Bernd Flemisch's avatar
Bernd Flemisch committed
383
  Compile the main file by typing \texttt{make tutorial\_coupled} and
384
  run the model as explained above.
Bernd Flemisch's avatar
Bernd Flemisch committed
385

386
  \item \textbf{Changing  the Shape of the Discrete Elements} \\
Thomas Fetzer's avatar
Thomas Fetzer committed
387
388
389
390
391
392
393
394
395
396
  In order to complete this exercise you need an external grid manager capable of handling
  simplex grids, like \texttt{ALUGrid} or \texttt{UGGrid}. If this is not the case,
  please skip this exercise.
  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/io/}, change for example from
  \texttt{cubegridcreator.hh} to \texttt{simplexgridcreator.hh}.
  For ALUGrid you have to change the ALUGrid type in line \ref{tutorial-coupled:set-grid-ALU} from \texttt{Dune::cube}
397
  to \texttt{Dune::simplex}.
398
  The shape of the employed elements can be visualized in ParaView by choosing \texttt{Surface with Edges}.
399

Bernd Flemisch's avatar
Bernd Flemisch committed
400
\item \textbf{Changing Fluids} \\
Thomas Fetzer's avatar
Thomas Fetzer committed
401
402
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:
403
\begin{enumerate}
Thomas Fetzer's avatar
Thomas Fetzer committed
404
405
406
407
408
 \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<Scalar>},
  as a second template argument after the data type \texttt{<Scalar>}, i.e.
  \texttt{Dumux::Brine<Scalar, Dumux::H2O<Scalar>>}. The file is located in the folder \texttt{dumux/material/components/}.
409
  Try to include the file and select the component as the wetting phase via the property system.
Thomas Fetzer's avatar
Thomas Fetzer committed
410
 \item DNAPL:
411
  Now let's include a DNAPL (\textbf{d}ense \textbf{n}on-\textbf{a}queous \textbf{p}hase \textbf{l}iquid)
Thomas Fetzer's avatar
Thomas Fetzer committed
412
413
  which 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.
414
\end{enumerate}
Thomas Fetzer's avatar
Thomas Fetzer committed
415
416
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
417
\texttt{/dumux/material/components} and read chapter~\ref{sec:fluidframework}.
Bernd Flemisch's avatar
Bernd Flemisch committed
418

419
\item \textbf{Use a Full-Fledged Fluid System} \\
Thomas Fetzer's avatar
Thomas Fetzer committed
420
421
422
423
\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} --
424
the same as to cancel the comment later on.\\
Thomas Fetzer's avatar
Thomas Fetzer committed
425
426
Now include the file \texttt{fluidsystems/h2oairfluidsystem.hh} in the material
folder, and set a type property \texttt{FluidSystem} (see line \ref{tutorial-coupled:set-fluidsystem})
427
428
429
430
431
with the appropriate type, which is either:\\
 \texttt{Dumux::FluidSystems::H2OAir<typename GET\_PROP\_TYPE(TypeTag, Scalar)>}\\
or in the \Dumux tongue\\
 \texttt{Dumux::H2OAirFluidSystem<TypeTag>}
\\
Thomas Fetzer's avatar
Thomas Fetzer committed
432
433
434
435
436
437
438
439
440
441
442
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();}.
Remember that the constructor function always has the same name as the respective
class, i.e. \texttt{TutorialProblemCoupled(..)}.\\
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.\\
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
443
444

\item \textbf{Changing Constitutive Relations} \\
Thomas Fetzer's avatar
Thomas Fetzer committed
445
446
  Use an unregularized linear law with an entry pressure of $p_e = \unit[0.0]{Pa}$
  and maximal capillary pressure of e.g. $p_{c_{max}} = \unit[2000.0]{Pa}$ instead of using a
447
 regularized Brooks-Corey law for the
Thomas Fetzer's avatar
Thomas Fetzer committed
448
449
450
451
452
453
  relative permeability and for the capillary pressure saturation relationship.
  To do that you have
  to change the material law property (line \ref{tutorial-coupled:eff2abs}) in
  \texttt{tutorialspatialparams\_coupled.hh}. Leave the type definition of \texttt{Scalar} and remove
 the type definition of \texttt{BrooksAndCorey} in the private section of
 the property definition. Exchange the \texttt{EffToAbsLaw} with the \texttt{LinearMaterial} law type in the
454
public section.
Thomas Fetzer's avatar
Thomas Fetzer committed
455
 You can find the material laws in the folder
Melanie Darcis's avatar
Melanie Darcis committed
456
  \verb+dumux/material/fluidmatrixinteractions+. The necessary parameters
457
458
of the linear law and the respective \texttt{set}-functions can be found
 in the file \\
459
460
 \verb+dumux/material/fluidmatrixinteractions/2p/linearmaterialparams.hh+.\\
Call the \texttt{set}-functions from the constructor of the \texttt{tutorialspatialparams\_coupled.hh}.
Thomas Fetzer's avatar
Thomas Fetzer committed
461

Bernd Flemisch's avatar
Bernd Flemisch committed
462
463
464
\item \textbf{Heterogeneities}  \\
  Set up a model domain with the soil properties given in Figure
  \ref{tutorial-coupled:exercise1_d}. Adjust the boundary conditions
465
  so that water is again flowing from the left to the right of the
466
\begin{figure}[ht]
Bernd Flemisch's avatar
Bernd Flemisch committed
467
\centering
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
\begin{tikzpicture}[>=latex]
  % basic sketch
  \fill [dumuxBlue](3,0) rectangle ++(3,3);
  \draw (3,0) -- ++(0,3);
  \draw (0,0) rectangle ++(6,3);
  % arrows
  \draw [|<->|](0,-0.3) -- ++(6,0);
  \node at (3,-0.3)[anchor=north]{$\unit[600]{m}$};
  \draw [|<->|](-0.3,0) -- ++(0,3);
  \node at(-0.3,1.5)[anchor=east]{$\unit[300]{m}$};
  % labels
  \node [anchor=west] at (0.2,1.5){$\mathbf{K}=\unit[10^{-8}]{m^2}$};
  \node [anchor=west] at (0.2,1){$\phi=0.15$};
  \node [anchor=west] at (3.2,1.5){$\mathbf{K}=\unit[10^{-9}]{m^2}$};
  \node [anchor=west] at (3.2,1){$\phi=0.3$};
Thomas Fetzer's avatar
Thomas Fetzer committed
483
484
485
\end{tikzpicture}
\caption{Exercise 1g: Set-up of a model domain with a heterogeneity. Grid
spacing: $\Delta x = \unit[20]{m}$ $\Delta y = \unit[20]{m}$.}\label{tutorial-coupled:exercise1_d}
Bernd Flemisch's avatar
Bernd Flemisch committed
486
\end{figure}
487
domain. You can use the fluids of exercise 1b.\\
Thomas Fetzer's avatar
Thomas Fetzer committed
488
489
\textbf{Hint:} The current position of the control volume can be obtained
using \texttt{element\allowbreak.geometry()\allowbreak.corner(scvIdx)}, which
490
returns a vector of the global coordinates of the current position.\\
491
When does the front cross the material border? In ParaView, the
492
animation view (\textit{View} $\rightarrow$ \textit{Animation
493
494
  View}) is a convenient way to get a rough feeling of the time-step
sizes.
Bernd Flemisch's avatar
Bernd Flemisch committed
495
496
497
\end{enumerate}

\subsubsection{Exercise 2}
498
For this exercise you should create a new problem file analogous to
499
the file \texttt{tutorialproblem\_coupled.hh} (e.g. with the name
500
501
\texttt{ex2\_tutorialproblem\_coupled.hh} and new spatial parameters \texttt{ex2\_tutorialspatialparams\_coupled.hh},
just like \texttt{tutorialspatialparams\_coupled.hh}.
502

503
504
505
506
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
507
\ref{tutorial-coupled:guardian2}
508
509
in the header files (e.g. change
\mbox{\texttt{DUMUX\_TUTORIALPROBLEM\_COUPLED\_HH}} to\\
Thomas Fetzer's avatar
Thomas Fetzer committed
510
\mbox{\texttt{DUMUX\_EX2\_TUTORIALPROBLEM\_COUPLED\_HH}}). Include the new problem file in \texttt{tutorial\_coupled.cc}.
511
Besides adjusting the guardian macros, the new problem file should define and
512
use a new type tag for the problem as well as a new problem class
Thomas Fetzer's avatar
Thomas Fetzer committed
513
514
515
e.g. \mbox{\texttt{Ex2TutorialProblemCoupled}}. The type tag definition has
to be adjusted in \texttt{tutorial\_coupled.cc} too (see line \ref{tutorial-coupled:set-type-tag}).
Similarly adjust your new spatial parameters file. If you are using Eclipse there is
516
a very helpful function called \texttt{Refactor} which you can use to change
Thomas Fetzer's avatar
Thomas Fetzer committed
517
518
all similar variables or types in your current file in one go. Just place the
cursor at the variable or type you want to change
519
and use the \texttt{Refactor} $\rightarrow$ \texttt{Rename} functionality. Make sure to assign your
520
newly defined spatial parameter class to the
521
\texttt{SpatialParams} property for the new
Thomas Fetzer's avatar
Thomas Fetzer committed
522
type tag.
523

524
525
526
527
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
528
529
saturated with water and the pressure is $p_w = \unit[5 \times
10^5]{Pa}$. Oil infiltrates from the left side. Create a grid
530
with $20$ cells in $x$-direction and $10$ cells in $y$-direction. The
531
532
simulation time should be set to $\unit[10^6]{s}$ with an
initial time-step size of $\unit[100]{s}$. Then, you can compile the program.
Bernd Flemisch's avatar
Bernd Flemisch committed
533
534


535
\begin{figure}[ht]
Bernd Flemisch's avatar
Bernd Flemisch committed
536
\centering
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
\begin{tikzpicture}[scale=0.7,>=latex]
  % basic sketch
  \draw (0,0) rectangle ++(10,5);
  \draw [fill=dumuxYellow] (2.5,1.5) rectangle ++(5,2);
  % arrows
  \draw[|<->|] (-0.2,0) -- ++(0,5);
  \node [anchor=east] at (-0.2,2.5){$\unit[50]{m}$};
  \draw[|<->|] (0,5.2) -- ++(10,0);
  \node [anchor=south] at (5,5.2){$\unit[100]{m}$};
  \draw[|<->|] (2.3,1.5) -- ++(0,2);
  \node [anchor=east] at (2.3,2.5){$\unit[20]{m}$};
  \draw[|<->|] (2.3,0) -- ++(0,1.5);
  \node [anchor=east] at (2.3,0.75){$\unit[15]{m}$};
  \draw[|<->|] (2.5,3.7) -- ++(5,0);
  \node [anchor=south] at (5,3.7){$\unit[50]{m}$};
  \draw[|<->|] (7.5,3.7) -- ++(2.5,0);
  \node [anchor=south] at (8.25,3.7){$\unit[25]{m}$};
  % labels
  \draw [dashed] (11,3) rectangle ++(7,3.5);
  \node [anchor=south west] at (11,5.4){$\mathbf{K} = \unit[10^{-7}]{m^2}$};
  \node [anchor=south west] at (11,4.6){$\phi = 0.2$};
  \node [anchor=south west] at (11,3.8){\textsc{Brooks-Corey Law}};
  \node [anchor=south west] at (11,3.0){$\lambda = 1.8, p_e = \unit[1000]{Pa}$};
  \draw [->] (11,4) -- (9.5,2.5);
  \draw [dashed] (11,-1) rectangle ++(7,3.5);
  \node [anchor=south west] at (11,1.4){$\mathbf{K} = \unit[10^{-9}]{m^2}$};
  \node [anchor=south west] at (11,0.6){$\phi = 0.15$};
  \node [anchor=south west] at (11,-0.2){\textsc{Brooks-Corey Law}};
  \node [anchor=south west] at (11,-1.0){$\lambda = 2, p_e = \unit[1500]{Pa}$};
  \draw [->] (11,1.5) -- (7,2.5);
\end{tikzpicture}
Bernd Flemisch's avatar
Bernd Flemisch committed
568
569
570
\caption{Set-up of the model domain and the soil parameters}\label{tutorial-coupled:ex2_Domain}
\end{figure}

571
\begin{figure}[ht]
Bernd Flemisch's avatar
Bernd Flemisch committed
572
\centering
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
\begin{tikzpicture}[scale=0.7,>=latex]
  % basic sketch
  \fill [pattern=north west lines] (0,0) rectangle ++(10,-0.25);
  \fill [pattern=north west lines] (0,5) rectangle ++(10,0.25);
  \draw (0,0) rectangle ++(10,5);
  \draw [fill=dumuxYellow] (2.5,1.5) rectangle ++(5,2);
  \foreach \y in {0.5,1.5,...,4.5}
    \draw [->](10.2,\y) -- ++(0.8,0);
  % labels
  \node [anchor=south] at (5,5.25){no flow};
  \node [anchor=north] at (5,-0.25){no flow};
  \node [anchor=west] at(11,2){$q_n = 0$};
  \node [anchor=west] at(11,3){$q_w = \unitfrac[2 \cdot 10^{-4}]{kg}{m^2 s}$};
  \node [anchor=west] at(-4,2){$S_n = 1$};
  \node [anchor=west] at(-4,3){$p_w = \unit[5 \cdot 10^5]{Pa}$};
\end{tikzpicture}
Bernd Flemisch's avatar
Bernd Flemisch committed
589
590
591
\caption{Boundary Conditions}\label{tutorial-coupled:ex2_BC}
\end{figure}

592
\begin{itemize}
593
 \item Increase the simulation time to e.g. $\unit[4\times 10^7]{s}$. Investigate the saturation: Is the value range reasonable?
594
595
596
 \item What happens if you increase the resolution of the grid?
\end{itemize}

597
\subsubsection{Exercise 3: Parameter File Input}
598
599
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
600
via \textit{parameter input files}.
601

602
603
In the code, parameters can be read via the macro
\texttt{GET\_RUNTIME\_PARAM(TypeTag, Scalar,
Thomas Fetzer's avatar
Thomas Fetzer committed
604
MyWonderfulGroup.MyWonderfulParameter);}. In this exercise we will explore the possibilities of the
605
parameter file. For this we take a look at the file \texttt{ex3\_tutorial\_coupled.input} in the \texttt{solutions\_coupled} folder.
Thomas Fetzer's avatar
Thomas Fetzer committed
606
607
608
609
610
611
612
613
Besides the parameters which you already used in the parameter file above,
there are parameters which can be used to control the
Newton and the Linear solver (groups: \texttt{Newton} and \texttt{LinearSolver}).
Run-time parameters used in the problem or spatial parameters classes
can also be set with the respective group names (\texttt{Problem} and \texttt{SpatialParams})
in the parameter file. For the latter parameters to be included in the program
they have to be assigned in the problem or spatial parameters constructor. This
can be done as shown in the files \texttt{ex3\_tutorialproblem\_coupled.diff}
614
615
and \texttt{ex3\_tutorialspatialparams\_coupled.diff} in the \texttt{solutions\_coupled} folder. Add some (for
example \texttt{Newton.MaxSteps} and \texttt{Problem.EnableGravity}) to the
616
parameter file \texttt{tutorial\_coupled.input} and observe what
617
happens if they are modified.
618

619
\subsubsection{Exercise 4: Create a New Component}
620
621
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
622
existing components in the directory \verb+/dumux/material/components+). \\
Bernd Flemisch's avatar
Bernd Flemisch committed
623
Use benzene as a new fluid and run the model of Exercise 2 with water
624
625
and benzene. Benzene has a density of $\unitfrac[889.51]{kg}{m^3}$ and
a viscosity of $\unit[0.00112]{Pa \, s}$.
626
627

\subsubsection{Exercise 5: Time Dependent Boundary Conditions}
Thomas Fetzer's avatar
Thomas Fetzer committed
628
629
630
In this exercise we want to investigate the influence of time dependent boundary
conditions. For this, redo the steps of exercise 2 and create a new problem and
spatial parameters file.
631
632
633
634

After this, change the run-time parameters so that they match the
domain described by figure \ref{tutorial-coupled:ex5_Domain}. Adapt
the problem class so that the boundary conditions are consistent with
Thomas Fetzer's avatar
Thomas Fetzer committed
635
636
637
638
639
640
figure \ref{tutorial-coupled:ex5_BC}. Here you can see the time dependence of the
nonwetting saturation, where water infiltrates only during $\unit[10^5]{s}$ and
$\unit[4 \cdot 10^5]{s}$. To implement these time dependencies you need the actual
time $t_{n+1}=t_n + \Delta t$ and the endtime of the simulation. For this you can
use the methods \texttt{this->timeManager().time()}, \texttt{this->timeManager().timeStepSize()}
and \texttt{this->timeManager().endTime()}.
641
642
643
644

Initially, the domain is fully saturated with oil and the pressure is $p_w = 2 \times
10^5\,\text{Pa}$.  Water infiltrates from the left side. Create a grid
with $100$ cells in $x$-direction and $10$ cells in $y$-direction. The
645
646
simulation time should be set to $\unit[5 \cdot 10^5]{s}$ with an
initial time-step size of $\unit[10]{s}$. To avoid too big time-step sizes
Thomas Fetzer's avatar
Thomas Fetzer committed
647
648
you should set the parameter \texttt{MaxTimeStepSize} for the group \texttt{TimeManager}
(in your input file) to $\unit[1000]{s}$.  Then, you can compile the program.
649
650
651
652
653
654
655
656
657
658
659
660
661
662

\begin{figure}[ht]
\centering
\begin{tikzpicture}[scale=0.7,>=latex]
  % basic sketch
  \fill [pattern=north west lines] (0,0) rectangle ++(10,-0.25);
  \fill [pattern=north west lines] (0,5) rectangle ++(10,0.25);
  \draw (0,0) rectangle ++(10,5);
  \foreach \y in {0.5,1.5,...,4.5}
    \draw [->](10.2,\y) -- ++(0.8,0);
  % arrows
  \draw[|<->|] (-0.2,0) -- ++(0,5);
  \node [anchor=east] at (-0.2,1.5){$\unit[50]{m}$};
  \draw[|<->|] (0,5.35) -- ++(10,0);
Thomas Fetzer's avatar
Thomas Fetzer committed
663
  \node [anchor=south] at (2.5,5.2){$\unit[100]{m}$};
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
  % labels
  \node [anchor=south] at (5,5.25){no flow};
  \node [anchor=north] at (5,-0.25){no flow};
  \node [anchor=west] at(11,2){$q_n = \unitfrac[1 \cdot 10^{-3}]{kg}{m^2 s}$};
  \node [anchor=west] at(11,3){$q_w = 0$};
  \node [anchor=west] at(-4,2){$S_n(t)$};
  \node [anchor=west] at(-4,3){$p_w = \unit[2 \cdot 10^5]{Pa}$};
  \node [anchor=south west] at (2.5,3.4){$\mathbf{K} = \unit[10^{-7}]{m^2}$};
  \node [anchor=south west] at (2.5,2.6){$\phi = 0.2$};
  \node [anchor=south west] at (2.5,1.8){\textsc{Brooks-Corey Law}};
  \node [anchor=south west] at (2.5,1.0){$\lambda = 2, \; p_e = \unit[500]{Pa}$};
\end{tikzpicture}
\caption{Set-up of the model domain and the soil parameters}\label{tutorial-coupled:ex5_Domain}
\end{figure}

\begin{figure}[ht]
\centering
\begin{tikzpicture}[scale=0.9,>=latex]
    % Draw axes
    \draw [<->,thick] (0,6) node (yaxis) [above] {$S_n$}
        |- (11,0) node (xaxis) [right] {time\,[s]};
    \draw plot[smooth,samples=100,domain=0:1] (6*\x + 2 ,{5*sin((\x+1)*pi r)+5});
	\draw [-] (0,5) -- (2,5);
	\draw [-] (8,5) -- (10,5);
    \draw [dashed] (2,5) -- (2,0);
    \draw [dashed] (8,5) -- (8,0);
    \draw [dashed] (10,5) -- (10,0);
Thomas Fetzer's avatar
Thomas Fetzer committed
691

692
693
694
695
696
697
698
699
700
    % axes labeling
    \draw [-] (-0.1,5) -- (0.1,5);
    \node [anchor=west] at(-0.5,5){$1$};
    \draw [-] (-0.1,0) -- (0.1,0);
    \node [anchor=west] at(-0.5,0){$0$};
    \draw [-] (2,0.1) -- (2,-0.1);
    \node [anchor=west] at(1.5,-0.4){$1\cdot10^{5}$};
    \draw [-] (8,0.1) -- (8,-0.1);
    \node [anchor=west] at(7.5,-0.4){$4\cdot10^{5}$};
Thomas Fetzer's avatar
Thomas Fetzer committed
701
702
703
    \draw [-] (10,0.1) -- (10,-0.1);
    \node [anchor=west] at(9.5,-0.4){$5\cdot10^{5}$};

704
    \node [anchor=base] at (5,3){$1-\sin\left(\pi\frac{\text{time}-10^5}{3\cdot 10^5 }\right)$};
705
706
707
708
709
\end{tikzpicture}
\caption{Time Dependent Boundary Conditions}\label{tutorial-coupled:ex5_BC}
\end{figure}

\begin{itemize}
Thomas Fetzer's avatar
Thomas Fetzer committed
710
711
 \item Open ParaView and plot the values of $S_n$ at time $\unit[5 \cdot 10^5]{s}$
 over the $x$-axis.\\ (\texttt{Filter->Data Analysis->Plot Over Line})
712
713
 \item What happens without any time-step restriction?
\end{itemize}
Bernd Flemisch's avatar
Bernd Flemisch committed
714

715
\clearpage \newpage
Thomas Fetzer's avatar
Thomas Fetzer committed
716
%%% Local Variables:
717
718
%%% mode: latex
%%% TeX-master: "dumux-handbook"
Thomas Fetzer's avatar
Thomas Fetzer committed
719
%%% End: