tutorial-decoupled.tex 20.2 KB
Newer Older
Bernd Flemisch's avatar
Bernd Flemisch committed
1
2
3
4
5
6
7
8
9
\section[Decoupled model]{Solving a problem using a Decoupled Model}\label{tutorial-decoupled}
The process of solving a problem using \Dumux can be roughly divided into four parts: 
\begin{enumerate}
 \item The geometry of the problem and correspondingly a grid have to be defined.
 \item Material properties and constitutive relationships have to be defined.
 \item Boundary conditions as well as initial conditions have to be defined.
 \item A suitable model has to be chosen.
\end{enumerate}

10
11
12
13
14
15
16
17
18
19
20
21
22
23
In contrast to the last section, we now apply a decoupled solution procedure, a
so-called \textit{IMPET} (\textit{IM}plicit \textit{P}ressure \textit{E}xplicit 
\textit{T}ransport) algorithm. This means that the pressure equation is first 
solved using an implicit method. The resulting velocities are then used to solve
a transport equation explicitly.\\
In this tutorial, pure fluid phases are solved with a finite volume discretization
of both pressure- and transport step. Primary variables, according to default
settings of the model, are the pressure and the saturation of the wetting phase.

The problem which is solved in this tutorial is illustrated in figure 
\ref{tutorial-decoupled:problemfigure}. A rectangular domain with now flow 
boundaries on the top and at the bottom, which is initially saturated with oil, 
is considered. Water infiltrates from the left side into the domain. Gravity 
effects are neglected.
Bernd Flemisch's avatar
Bernd Flemisch committed
24
25
26
27
28
29
30
31
32
33
34
35

\begin{figure}[h]
\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_w = 1$}
\psfrag{S_n_initial = 0}{\textcolor{white}{$\mathbf{S_{w_{initial}} = 0}$}}
\psfrag{q_w = 0 [kg/m^2s]}{$q_w = 0$ $\left[\frac{\textnormal{kg}}{\textnormal{m}^2 \textnormal{s}}\right]$}
36
\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
37
38
39
40
41
\centering
\includegraphics[width=0.9\linewidth,keepaspectratio]{EPS/tutorial-problemconfiguration}
\caption{Geometry of the tutorial problem with initial and boundary conditions.}\label{tutorial-decoupled:problemfigure}
\end{figure}

42
Listing \ref{tutorial-deoucpled:mainfile} shows how the main file, which has to be executed, has to be set up, if the problem described above is to be solved using a decoupled model. This main file can be found in the directory \texttt{/tutorial} of the stable part of \Dumux.
Bernd Flemisch's avatar
Bernd Flemisch committed
43
44
45

\begin{lst}[File tutorial/tutorial\_decoupled.cc]\label{tutorial-deoucpled:mainfile} \mbox{}
\lstinputlisting[basicstyle=\ttfamily\scriptsize,numbers=left, 
46
numberstyle=\tiny, numbersep=5pt, firstline=28]{../../tutorial/tutorial_decoupled.cc}
Bernd Flemisch's avatar
Bernd Flemisch committed
47
48
\end{lst}

49
First, from line \ref{tutorial-decoupled:include-begin} to line
50
\ref{tutorial-decoupled:include-end} the \Dune and \Dumux files containing
51
52
53
54
55
56
57
58
essential functions and classes are included.

At line \ref{tutorial-decoupled:set-type-tag} the type tag of the
problem which is going to be simulated is set. All other data types
can be retrieved by the \Dumux property system and only depend on this
single type tag. Retrieving them is done between line
\ref{tutorial-decoupled:retrieve-types-begin} and
\ref{tutorial-decoupled:retrieve-types-end}. For an introduction to the
59
property system, see section \ref{sec:propertysystem}.
60
61

The first thing which should be done at run time is to initialize the
62
message passing interface using \Dune's \texttt{MPIHelper} class. Line
63
\ref{tutorial-decoupled:init-mpi} is essential if the simulation is
64
65
66
intended to be run on more than one processor at the same time. Next,
the command line arguments are parsed starting at line
\ref{tutorial-decoupled:parse-args-begin} until line
67
\ref{tutorial-decoupled:parse-args-end}. In this case, we check weather and
68
69
70
71
72
73
74
at which time a previous run of the simulation should be restarted, and we
parse the time when the simulation ends. As the maximum time-step in the 
sequential model is strictly bound by a CFL-criterion, the first time-step
size is initialized with the simulation time.

After this, a grid is created on line \ref{tutorial-decoupled:create-grid} 
and the problem is instantiated with information about the grid
75
(via its leaf grid view) on line \ref{tutorial-decoupled:instantiate-problem}.
76
If demanded, on line \ref{tutorial-decoupled:mainRestart} a state written to
77
78
disk by a previous simulation run is restored on request by the user.
Finally, the time manager controlling the simulation run is instantiated 
79
with the start parameters in line \ref{tutorial-decoupled:initTimeManager}
80
and the simulation proceedure is started by the time manager at line
81
\ref{tutorial-decoupled:execute}.
82
83
84
85
86



\subsection{The problem class} \label{decoupled_problem}

87
88
89
90
91
92
93
When solving a problem using \Dumux, the most important file is the
so-called \textit{problem file} as shown in listing
\ref{tutorial-decoupled:problemfile} of
\texttt{tutorialproblem\_decoupled.hh}.

\begin{lst}[File tutorial/tutorialproblem\_decoupled.hh]\label{tutorial-decoupled:problemfile} \mbox{}
\lstinputlisting[basicstyle=\ttfamily\scriptsize,numbers=left,
94
numberstyle=\tiny, numbersep=5pt, firstline=28]{../../tutorial/tutorialproblem_decoupled.hh}
95
96
97
98
99
100
101
102
\end{lst}

First, both \Dune  grid handlers and the decoupled model of \Dumux 
have to be included. Then, a new type tag is created for the problem 
on line \ref{tutorial-decoupled:create-type-tag}.  In this case, the 
new type tag inherits all properties defined for the \texttt{DecoupledTwoP} 
type tag, which means that for this problem the two-phase decoupled approach
is chosen as discretization scheme (defined via the include in line 
103
\ref{tutorial-decoupled:parent-problem}). On line \ref{tutorial-decoupled:set-problem}, 
104
a problem class is attached to the new type tag, while the grid which
105
is going to be used is defined on line \ref{tutorial-decoupled:set-grid-type} --
106
in this case an \texttt{SGrid} is created.  Since in \Dune, there's no uniform
107
108
109
mechanism to allocate grids, the \texttt{Grid} property also contains
a static \texttt{create()} method which provides just that: From line 
\ref{tutorial-decoupled:grid-begin} to \ref{tutorial-decoupled:grid-end}, 
110
111
112
the geometry is defined and the grid is generated. The three variables of 
Type \texttt{Dune::FieldVector} define the lower left corner of the domain 
(\texttt{L}), the upper right corner of the domain (\texttt{H}) and the number 
113
of cells in $x$ and $y$ direction (\texttt{N}). The grid of type 
114
\texttt{Dune::SGrid} is then generated in line \ref{tutorial-decoupled:grid-end}. 
115
For more information about the \Dune grid interface, the different grid types 
116
that are supported and the generation of different grids it is referred to 
117
118
the \textit{Dune Grid Interface HOWTO} \cite{DUNE-HP}. 

119
120
121
Next, we select the material of the simulation: In case of a pure two-phase
model, each phase is a bulk fluid, and the complex (compositional) fluidsystems
do not need to be used. However, they can be used (see exercise 1 \ref{dec-ex1-fluidsystem}). 
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
Instead, we use a simplified fluidsystem container that provides classes 
for liquid and gas phases, line \ref{tutorial-decoupled:2p-system-start} to 
\ref{tutorial-decoupled:2p-system-end}. These are linked to the appropriate 
chemical species in line \ref{tutorial-decoupled:wettingPhase} and 
\ref{tutorial-decoupled:nonwettingPhase}. For all parameters that depend 
on space, such as the properties of the soil, the specific spatial parameters 
for the problem of interest are specified in line
\ref{tutorial-decoupled:set-spatialparameters}. 

Now we arrive at some model parameters of the applied two-phase decoupled 
model. Line \ref{tutorial-decoupled:velocityFormulation} defines that the 
wetting phase velocity rather than e.g. a total velocity is used for the 
transport system. As we regard capillary pressure, a capillary diffusive 
term is regarded, selected in line \ref{tutorial-decoupled:DiffusivePart}.
Line \ref{tutorial-decoupled:cfl} assigns the CFL-factor to be used in the
137
simulation run. The final property on line \ref{tutorial-decoupled:gravity} 
138
139
140
141
142
143
is optional and tells the model not to use gravity.

After all necessary information is written into the property system and 
its namespace is closed in line \ref{tutorial-decoupled:propertysystem-end},
the problem class is defined in line \ref{tutorial-decoupled:def-problem}. 
As its property, the problem class itsself is also derived from a parent, 
144
\texttt{IMPESProblem2P}. The class constructor (line 
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
\ref{tutorial-decoupled:constructor-problem}) is able to hold two vectors,
which is not needed in this tutorial.

Besides the definition of the boundary and initial conditions (discussed in 
subsection \label{decoupled-problem:boundary}), the problem class also contains
general information about the current simulation. First, the name used by
the \texttt{VTK-writer} to generate output is defined in the method of line
\ref{tutorial-decoupled:name}, and line \ref{tutorial-decoupled:restart} indicates
weather restart files are written. As decoupled schemes usually feature small 
timesteps, the method controlling the output in line \ref{tutorial-decoupled:output}
is very useful. The divisor of the modulo operation defines after how many timesteps
output should be written out -- the default ``1'' resembles output after each 
step.

The following methods all have in common that they may be dependent on space.
Hence, they all feature a common argument list:
\begin{itemize}
 \item \texttt{globalPos}: A vector holding the global Coordinates.
 \item \texttt{element} or \texttt{intersection}: Input for an iterator, that is 
    depending weather the parameter of the method is defined in an element, such as 
    initial values, or on an intersection, such as a boundary condition.
\end{itemize}
In the following, there are the methods for general parameters, source- or
sinkterms, boundary conditions (lines \ref{tutorial-decoupled:bctypePress} to
169
\ref{tutorial-decoupled:neumann}) and initial values for the transported
170
171
quantity in line \label{tutorial-decoupled:initSat}. For more information
on the functions, it is referred to the documentation in the code.
Bernd Flemisch's avatar
Bernd Flemisch committed
172

173
\subsection{The definition of the parameters that are dependent on space}\label{tutorial-decoupled:description-spatialParameters}
Bernd Flemisch's avatar
Bernd Flemisch committed
174

175
176
Listing \ref{tutorial-decoupled:spatialparametersfile} shows the file
\verb+tutorialspatialparameters_decoupled.hh+:
Bernd Flemisch's avatar
Bernd Flemisch committed
177

178
179
\begin{lst}[File tutorial/tutorialspatialparameters\_decoupled.hh]\label{tutorial-decoupled:spatialparametersfile} \mbox{}
\lstinputlisting[basicstyle=\ttfamily\scriptsize,numbers=left,
180
numberstyle=\tiny, numbersep=5pt, firstline=26]{../../tutorial/tutorialspatialparameters_decoupled.hh}
Bernd Flemisch's avatar
Bernd Flemisch committed
181
\end{lst}
182
183
184
185
186
As this file only slightly differs from the coupled version, it is referred to 
chapter \ref{tutorial-coupled:description-spatialParameters} for explanations.
However, as a standard Finite-Volume--scheme is used, in contrast to the box-method
in the coupled case, the argument list here is the same as for the problem 
functions:
Bernd Flemisch's avatar
Bernd Flemisch committed
187
\begin{itemize}
188
189
190
 \item \texttt{globalPos}: A vector holding the global Coordinates.
 \item \texttt{element}: Input for an element iterator, providing access
	to the current element of interest.
Bernd Flemisch's avatar
Bernd Flemisch committed
191
192
193
194
\end{itemize}

\subsection{Exercise}
\label{tutorial-deoucpled:exercises}
195
196
197
The following exercises will give you the opportunity to learn how you can change 
soil parameters, boundary conditions and fluid properties in \Dumux and to play along 
with the decoupled modelling framework.
Bernd Flemisch's avatar
Bernd Flemisch committed
198
199
200
201
202

\subsubsection{Exercise 1}
\renewcommand{\labelenumi}{\alph{enumi})}
For Exercise 1 you only have to make some small changes in the tutorial files.
\begin{enumerate}
203
\item \textbf{Altering output}
204
To get an impression what the results should look like you can first run the original version of the decoupled tutorial model by typing  \texttt{./tutorial\_decoupled 1e5}. The number behind the simulation name defines the timespan of the simulation run in seconds. For the visualisation with paraview please refer to \ref{quick-start-guide}.\\
205
As you can see, the simulation creates roughly 150 output files. To reduce these to perform longer simulations, change the method responsible for output (line \ref{tutorial-decoupled:output} in the file \texttt{tutorialproblem\_decoupled}) to write an output only every 20 timesteps by changeing the divisor. Compile the main file by typing \texttt{make tutorial\_decoupled} and run the model. Now, run the simulation for 5e5 seconds.
206

Bernd Flemisch's avatar
Bernd Flemisch committed
207
208
\item \textbf{Changing the Model Domain and the Boundary Conditions} \\
Change the size of the model domain so that you get a rectangle
209
with edge lengths of x = 300 m \\  and y = 300 m and with discretisation lengths of  $\Delta \text{x} = 20$ m and $\Delta \text{y} = 10$ m. \\
210
Change the boundary conditions in the file \texttt{tutorialproblem\_decoupled.hh} so that water enters from the bottom and oil flows out at the top boundary. The right and the left boundary should be closed for water and oil fluxes.  \\
Bernd Flemisch's avatar
Bernd Flemisch committed
211
212

\item \textbf{Changing Fluids} \\
213
214
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:
\begin{enumerate}
215
 \item Brine: The class \texttt{Dumux::Brine} acts as a adapter to the fluid system that alters a pure water class by adding some salt. 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>} as a template argument (be aware to use the complete water class with its own template parameter).
216
 \item DNAPL: A standard set of chemical substances, such as Oil and Brinde, is already included (via a list of \texttt{\#include ..} commandos) and hence easy accessible by default. This is not the case for the class \texttt{Dumux::SimpleDNAPL}, however, which is located in the folder \texttt{dumux/material/components/}. Try to include the file as well as select the component via the property system.
217
218
219
\end{enumerate}
If you want to take a closer look how the fluid classes are defined and which substances are already available please browse through the files in the directory
\texttt{/dumux/material/components}.
Bernd Flemisch's avatar
Bernd Flemisch committed
220

221
222
\item \textbf{Use the \Dumux fluid system}\label{dec-ex1-fluidsystem} \\
As you have experienced in the coupled tutorial (chapter \ref{tutorial-decoupled}), \Dumux usually organises fluid mixtures via a \texttt{fluidsystem}. This is also possible for the decoupled models: Uncomment, as we want to reuse it later on, the lines \ref{tutorial-decoupled:2p-system-start} to \ref{tutorial-decoupled:2p-system-end} in the problem file. If you use eclipse, this can easily be done by pressing \textit{str + shift + 7}, the same shortcut works to cancel the comment later on.\\
223
Now include the file \texttt{fluidsystems/h2o\_n2\_system.hh} in the material folder, and set a property \texttt{FluidSystem} with the appropriate type, \texttt{Dumux::H2O\_N2\_System<TypeTag>}. However, the complicated fluidsystem uses tabularized fluid data, which need to be initialized in the constructor body of the current problem by adding \texttt{GET\_PROP\_TYPE(TypeTag, PTAG(FluidSystem))::init();}, hence using the initialization function of the applied fluidsystem. As an alternative, use a simpler version of water, e.g. \texttt{Dumux::SimpleH2O}, and apply it for the property \texttt{Components} with type \texttt{H2O}. The density of the gas is magnitudes smaller than that of oil, so please decrease the injection rate to $q_n = -3 \times 10^-4$ $\left[\frac{\textnormal{kg}}{\textnormal{m}^2 \textnormal{s}}\right]$. Also reduce the simultation duration to 2e4 seconds.\\
224
Please reverse the changes of this example, as we still use bulk phases and hence do not need such an extensive fluid system.
Bernd Flemisch's avatar
Bernd Flemisch committed
225
226
 
\item \textbf{Heterogeneities}  \\
227
Set up a model domain with the soil properties given in Figure \ref{tutorial-deoucpled:exercise1_d}. Adjust the boundary conditions so that water is again flowing from left to right.
Bernd Flemisch's avatar
Bernd Flemisch committed
228
229
230
231
232
233
234
235
236
237
238
\begin{figure}[h]
\psfrag{K1 =}{K $= 10^{-8}\text{ m}^2$}
\psfrag{phi1 =}{$\phi = 0.15$}
\psfrag{K2 =}{\textcolor{white}{K $= 10^{-9}\text{ m}^2$}}
\psfrag{phi2 =}{\textcolor{white}{$\phi = 0.3$}}
\psfrag{600 m}{600 m}
\psfrag{300 m}{300 m}
\centering
\includegraphics[width=0.5\linewidth,keepaspectratio]{EPS/exercise1_c.eps}
\caption{Exercise 1d: Set-up of a model domain a heterogeneity. $\Delta \text{x} = 20$ m $\Delta \text{y} = 20$ m.}\label{tutorial-deoucpled:exercise1_d}
\end{figure}
239
When does the front cross the material border? In paraview, the option \textit{View} $\rightarrow$ \textit{Animation View} is nice to get a rough feeling of the timestep sizes.
Bernd Flemisch's avatar
Bernd Flemisch committed
240
241
242
\end{enumerate}

\subsubsection{Exercise 2}
243
For this exercise you should create a new proplem file analogous to
244
245
the file \texttt{tutorialproblem\_decoupled.hh} (e.g. with the name 
\texttt{ex2\_tutorialproblem\_decoupled.hh} and new spatial parameters 
246
247
just like \texttt{tutorialspatialparameters\_decoupled.hh}. These files need to
be included in the file \texttt{tutorial\_decoupled.cc}. 
248

249
250
251
252
253
254
255
The new files should contain the definition of a new classes with 
names that relate to the file name, such as \texttt{Ex2TutorialProblemDecoupled}. 
Make sure that you also adjust the guardian
macros in lines \ref{tutorial-decoupled:guardian1} and \ref{tutorial-decoupled:guardian1}
 in the header files (e.g. change \\
\texttt{DUMUX\_TUTORIALPROBLEM\_DECOUPLED\_HH} to
\texttt{DUMUX\_EX2\_TUTORIALPROBLEM\_DECOUPLED\_HH}).  Besides also adjusting the guardian macros, 
256
the new problem file should define and use a new type tag for the problem as well as a new problem class
257
e.g. \texttt{Ex2TutorialProblemDecoupled}. Make sure you assign your newly defined spatial 
258
259
260
261
262
parameter class to the \texttt{SpatialParameters} property for the new 
type tag. 

After this, change the \texttt{create()} method of the \texttt{Grid}
property so that it matches the domain described
263
by figure \ref{tutorial-decoupled:ex2_Domain}. Adapt the problem class
264
so that the boundary conditions are consistent with figure
265
\ref{tutorial-decoupled:ex2_BC}. Initially the domain is fully saturated
266
with water and the pressure is $p_w = 2 \times 10^5 \text{Pa}$ . Oil
267
268
infiltrates from the left side. Create a grid with $20$ cells in
$x$-direction and $10$ cells in $y$-direction. The simulation time
269
should be set to $2e4 \text{s}$.
270
271
272
273
274

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.

Bernd Flemisch's avatar
Bernd Flemisch committed
275
276
277
278

\begin{figure}[h]
\psfrag{K1}{K $= 10^{-7}\text{ m}^2$}
\psfrag{phi1}{$\phi = 0.2$}
279
\psfrag{Lin}{Brooks Corey Law} 
280
\psfrag{Lin2}{$\lambda = 1.8$, $p_b = 100$}
Bernd Flemisch's avatar
Bernd Flemisch committed
281
282
283
\psfrag{K2}{K $= 10^{-9}\text{ m}^2$}
\psfrag{phi2}{$\phi = 0.15$}
\psfrag{BC1}{Brooks Corey Law} 
284
\psfrag{BC2}{$\lambda = 2$, $p_b = 500$}
Bernd Flemisch's avatar
Bernd Flemisch committed
285
286
287
288
289
290
291
292
\psfrag{H1y}{50 m}
\psfrag{H2y}{15 m}
\psfrag{H3y}{20 m}
\psfrag{L1x}{100 m}
\psfrag{L2x}{50 m}
\psfrag{L3x}{25 m}
\centering
\includegraphics[width=0.8\linewidth,keepaspectratio]{EPS/Ex2_Domain.eps}
293
\caption{Set-up of the model domain and the soil parameters}\label{tutorial-decoupled:ex2_Domain}
Bernd Flemisch's avatar
Bernd Flemisch committed
294
295
296
\end{figure}

\begin{figure}[h]
297
298
299
\psfrag{pw}{$p_w = 2 \times 10^5$ [\text{Pa}]}
\psfrag{S}{$S_w = 0.0$}
\psfrag{qw}{$q_w = 3 \times 10^{-4}$ [kg/$\text{m}^2$s]}
Bernd Flemisch's avatar
Bernd Flemisch committed
300
301
302
303
\psfrag{qo}{$q_n = 0.0$ [kg/$\text{m}^2$s]}
\psfrag{no flow}{no flow}
\centering
\includegraphics[width=0.8\linewidth,keepaspectratio]{EPS/Ex2_Boundary.eps}
304
\caption{Boundary Conditions}\label{tutorial-decoupled:ex2_BC}
Bernd Flemisch's avatar
Bernd Flemisch committed
305
306
\end{figure}

307
\begin{itemize}
308
309
310
 \item What happens if you increase the resolution of the grid? Hint: Paraview can visualize the timesteps via the ``Animation View'' (to be enabled unter the button \textit{View}).
 \item Set the CFL-factor to 1 and investigate the saturation: Is the value range reasonable?
 \item Further increase the CFL-factor to 2 and investigate the saturation.
311
312
313
\end{itemize}


Bernd Flemisch's avatar
Bernd Flemisch committed
314
\subsubsection{Exercise 3}
315
316
Create a new file for benzene called \texttt{benzene.hh} and implement
a new fluid system. (You may get a hint by looking at existing fluid 
317
systems in the directory \verb+/dumux/material/fluidsystems+.)
318
319
320
321

Use benzene as a new fluid and run the model 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}$. 
322
323
324
325
326

%%% Local Variables: 
%%% mode: latex
%%% TeX-master: "dumux-handbook"
%%% End: