tutorial-coupled.tex 27.8 KB
 Andreas Lauser committed Feb 13, 2012 1 \section[Fully-Implicit Model]{Solving a Problem Using a Fully-Coupled Model}\label{tutorial-coupled}  Bernd Flemisch committed Jul 16, 2010 2   Klaus Mosthaf committed Feb 09, 2012 3 The process of setting up a problem using \Dumux can be roughly divided into four parts:  Bernd Flemisch committed Jul 16, 2010 4 \begin{enumerate}  Klaus Mosthaf committed Feb 09, 2012 5  \item A suitable model has to be chosen.  Bernd Flemisch committed Jul 16, 2010 6  \item The geometry of the problem and correspondingly a grid have to be defined.  Klaus Mosthaf committed Oct 19, 2010 7  \item Material properties and constitutive relationships have to be selected.  Philipp Nuske committed Feb 14, 2012 8  \item Boundary conditions and initial conditions have to be specified.  Bernd Flemisch committed Jul 16, 2010 9 10 \end{enumerate}  Klaus Mosthaf committed Feb 09, 2012 11 The problem being solved in this tutorial is illustrated in Figure \ref{tutorial-coupled:problemfigure}.  Klaus Mosthaf committed Feb 23, 2012 12 A rectangular domain with no-flow boundaries on the top and on the bottom, which is initially saturated with oil, is considered.  Klaus Mosthaf committed Feb 09, 2012 13 Water infiltrates from the left side into the domain and replaces the oil. Gravity effects are neglected here.  Bernd Flemisch committed Jul 16, 2010 14   15 \begin{figure}[ht]  Bernd Flemisch committed Jul 16, 2010 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]$}  Melanie Darcis committed Dec 07, 2010 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 committed Jul 16, 2010 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}  Klaus Mosthaf committed Feb 09, 2012 32 The solved equations are the mass balances of water and oil:  Bernd Flemisch committed Jul 16, 2010 33 34 35 \begin{align} \label{massbalancewater} \frac {\partial (\phi \, S_{w}\, \varrho_{w})}{\partial t}  Philipp Nuske committed Dec 10, 2010 36  -  Melanie Darcis committed Nov 11, 2010 37  \nabla \cdot \left( \varrho_{w} \, \frac{k_{rw}}{\mu_{w}} \, \mathbf{K}\;\nabla p_w \right)  Bernd Flemisch committed Jul 16, 2010 38  -  Melanie Darcis committed Nov 11, 2010 39  q_w  Bernd Flemisch committed Jul 16, 2010 40 41 42 43  & = 0 \\ \label{massbalanceoil} \frac {\partial (\phi \, S_{o}\, \varrho_{o})}{\partial t}  Philipp Nuske committed Dec 10, 2010 44  -  Melanie Darcis committed Nov 11, 2010 45  \nabla \cdot \left( \varrho_{o} \, \frac{k_{ro}}{\mu_{o}} \, \mathbf{K}\;\nabla p_o \right)  Bernd Flemisch committed Jul 16, 2010 46  -  Melanie Darcis committed Nov 11, 2010 47  q_o  Bernd Flemisch committed Jul 16, 2010 48 49 50 51  & = 0 \end{align}  Andreas Lauser committed Feb 13, 2012 52 \subsection{The Main File}  Bernd Flemisch committed Jul 16, 2010 53   Klaus Mosthaf committed Feb 09, 2012 54 Listing \ref{tutorial-coupled:mainfile} shows the main application file  Klaus Mosthaf committed Oct 19, 2010 55 \texttt{tutorial/tutorial\_coupled.cc} for the coupled two-phase  Klaus Mosthaf committed Feb 09, 2012 56 model. This file has to be compiled and executed in order to solve the problem described  Bernd Flemisch committed Jul 16, 2010 57 58 59 above. \begin{lst}[File tutorial/tutorial\_coupled.cc]\label{tutorial-coupled:mainfile} \mbox{}  Christoph Grueninger committed Feb 24, 2012 60  \lstinputlisting[style=DumuxCode, numbersep=5pt, firstline=29]{../../tutorial/tutorial_coupled.cc}  Bernd Flemisch committed Jul 16, 2010 61 62 63 \end{lst} From line \ref{tutorial-coupled:include-begin} to line  Klaus Mosthaf committed Feb 23, 2012 64 \ref{tutorial-coupled:include-end} the required headers are included.  Bernd Flemisch committed Jul 16, 2010 65 66  At line \ref{tutorial-coupled:set-type-tag} the type tag of the  Andreas Lauser committed Feb 13, 2012 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  Klaus Mosthaf committed Feb 23, 2012 69 on this single type tag. For a more thorough introduction to the  Andreas Lauser committed Feb 13, 2012 70 71 \Dumux property system, see chapter~\ref{sec:propertysystem}.  Klaus Mosthaf committed Feb 23, 2012 72 After this, the default startup routine \texttt{Dumux::start()} is  Andreas Lauser committed Feb 14, 2012 73 74 called on line \ref{tutorial-coupled:call-start}. This function deals with parsing the command line arguments, reading the parameter file,  Klaus Mosthaf committed Feb 23, 2012 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  Andreas Lauser committed Feb 14, 2012 80 81 (\texttt{-ParameterName ParameterValue}), in the file specified by the \texttt{-parameterFile} argument, or if the latter is not specified,  Klaus Mosthaf committed Feb 23, 2012 82 83 in the file \mbox{\texttt{tutorial\_coupled.input}}. If a parameter is  Andreas Lauser committed Feb 14, 2012 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 committed Feb 24, 2012 90 \lstinputlisting[style=DumuxCode]{../../tutorial/tutorial_coupled.input}  Andreas Lauser committed Feb 14, 2012 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  Klaus Mosthaf committed Feb 23, 2012 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}.  Philipp Nuske committed Feb 14, 2012 100   Andreas Lauser committed Feb 13, 2012 101 102  \subsection{The Problem Class}  Bernd Flemisch committed Jul 16, 2010 103 104  When solving a problem using \Dumux, the most important file is the  Andreas Lauser committed Feb 13, 2012 105 106 so-called \textit{problem file} as shown in listing~\ref{tutorial-coupled:problemfile}.  Bernd Flemisch committed Jul 16, 2010 107 108  \begin{lst}[File tutorial/tutorialproblem\_coupled.hh]\label{tutorial-coupled:problemfile} \mbox{}  Christoph Grueninger committed Feb 24, 2012 109 \lstinputlisting[style=DumuxCode, numbersep=5pt, firstline=28]{../../tutorial/tutorialproblem_coupled.hh}  Bernd Flemisch committed Jul 16, 2010 110 111 \end{lst}  Klaus Mosthaf committed Oct 18, 2010 112 First, a new type tag is created for the problem in line  Bernd Flemisch committed Jul 16, 2010 113 \ref{tutorial-coupled:create-type-tag}. In this case, the new type  Andreas Lauser committed Feb 13, 2012 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  Philipp Nuske committed Feb 14, 2012 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  Andreas Lauser committed Feb 13, 2012 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  Klaus Mosthaf committed Feb 23, 2012 131 the command-line or in a parameter file.  Andreas Lauser committed Feb 13, 2012 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  Klaus Mosthaf committed Feb 23, 2012 136 assumes immiscibility of the phases, but requires the components  Andreas Lauser committed Feb 13, 2012 137 138 used for the wetting and non-wetting phases to be explicitly set. In this case, liquid water which uses the relations from  Andreas Lauser committed Nov 10, 2010 139 IAPWS'97~\cite{IAPWS1997} is chosen as the wetting phase on line  Andreas Lauser committed Feb 13, 2012 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 committed Jul 16, 2010 152 153 154  The problem class always has at least five methods: \begin{itemize}  Klaus Mosthaf committed Feb 09, 2012 155 156 \item A method \texttt{boundaryTypes()} specifying the type of boundary conditions at each vertex.  Bernd Flemisch committed Jul 16, 2010 157 \item A method \texttt{dirichlet()} specifying the actual values for  Andreas Lauser committed Feb 13, 2012 158  the \textsc{Dirichlet} conditions at each \textsc{Dirichlet} vertex.  Bernd Flemisch committed Jul 16, 2010 159 \item A method \texttt{neumann()} specifying the actual values for  Andreas Lauser committed Feb 13, 2012 160 161  the \textsc{Neumann} conditions, which are usually evaluated at the integration points of the \textsc{Neumann} boundary faces.  Klaus Mosthaf committed Feb 09, 2012 162 \item A method for source or sink terms called \texttt{source()}, usually evaluated at  Andreas Lauser committed Feb 13, 2012 163  the center of a control volume.  Bernd Flemisch committed Jul 16, 2010 164 \item A method called \texttt{initial()} for specifying the initial  Klaus Mosthaf committed Feb 09, 2012 165  conditions at each vertex.  Bernd Flemisch committed Jul 16, 2010 166 167 \end{itemize}  Andreas Lauser committed Feb 13, 2012 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:  Klaus Mosthaf committed Feb 25, 2011 171 172 \begin{description} \item [values:] A vector which stores the result of the method. What  Viktor Szukitsch committed Sep 01, 2011 173  the values in this vector mean is dependent on the method: For  Klaus Mosthaf committed Feb 25, 2011 174 175  \texttt{dirichlet()} it contains the actual values of the primary variables, for \texttt{boundaryTypes()} it contains the boundary  Klaus Mosthaf committed Feb 09, 2012 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  Andreas Lauser committed Feb 13, 2012 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  Klaus Mosthaf committed Feb 23, 2012 184  volumes. Be aware that the second parameter is a Dune grid entity  Andreas Lauser committed Feb 13, 2012 185  with the codimension \texttt{dim}.  Klaus Mosthaf committed Feb 25, 2011 186 187 \end{description}  Andreas Lauser committed Feb 13, 2012 188 189 To ensure that no boundaries are undefined, a small safeguard value \texttt{eps\_} is usually added when comparing spatial  Klaus Mosthaf committed Feb 23, 2012 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  Klaus Mosthaf committed Feb 25, 2011 192 193 smaller than a very small value \texttt{eps\_}.  Andreas Lauser committed Feb 13, 2012 194 Methods which make statements about boundary segments of the grid  Klaus Mosthaf committed Feb 23, 2012 195 (such as \texttt{neumann()}) are called with six arguments:  Bernd Flemisch committed Jul 16, 2010 196 \begin{description}  Klaus Mosthaf committed Feb 25, 2011 197 198 \item[values:] A vector \texttt{neumann()}, in which the mass fluxes per area unit over the boundary segment are specified.  Bernd Flemisch committed Jul 16, 2010 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.  Klaus Mosthaf committed Feb 09, 2012 203 \item[isIt:] The \texttt{Intersection} of the boundary segment as given by the grid.  Bernd Flemisch committed Jul 16, 2010 204 \item[scvIdx:] The index of the sub-control volume in  Klaus Mosthaf committed Feb 23, 2012 205  \texttt{fvElementGeometry} which is assigned to the boundary segment.  Bernd Flemisch committed Jul 16, 2010 206 207 208 209 \item[boundaryFaceIdx:] The index of the boundary face in \texttt{fvElementGeometry} which represents the boundary segment. \end{description}  Melanie Darcis committed Oct 13, 2010 210 Similarly, the \texttt{initial()} and \texttt{source()} methods  Andreas Lauser committed Feb 13, 2012 211 specify properties of control volumes and thus only get  Bernd Flemisch committed Jul 16, 2010 212 \texttt{values}, \texttt{element}, \texttt{fvElemGeom} and  Andreas Lauser committed Feb 13, 2012 213 \texttt{scvIdx} as arguments.  Bernd Flemisch committed Jul 16, 2010 214 215  In addition to these five methods, there might be some model-specific  Andreas Lauser committed Feb 13, 2012 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  Klaus Mosthaf committed Feb 23, 2012 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  Andreas Lauser committed Feb 13, 2012 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}. \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 committed Jul 16, 2010 242   Viktor Szukitsch committed Sep 01, 2011 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.  Klaus Mosthaf committed Oct 19, 2010 245   Klaus Mosthaf committed Feb 23, 2012 246 \subsection{Defining Spatially Dependent Parameters}\label{tutorial-coupled:description-spatialParameters}  Bernd Flemisch committed Jul 16, 2010 247   Andreas Lauser committed Feb 13, 2012 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  Klaus Mosthaf committed Feb 23, 2012 253 parameters are defined using a so-called \textit{spatial parameters}  Andreas Lauser committed Feb 13, 2012 254 class.  Bernd Flemisch committed Jul 16, 2010 255   Klaus Mosthaf committed Feb 23, 2012 256 If the box discretization is to used, the spatial parameters class  Andreas Lauser committed Feb 13, 2012 257 should be derived from the base class  258 \texttt{Dumux::BoxSpatialParams}. Listing  Klaus Mosthaf committed Feb 23, 2012 259 \ref{tutorial-coupled:spatialparametersfile} shows the file \\  260 \verb+tutorialspatialparams_coupled.hh+:  Bernd Flemisch committed Jul 16, 2010 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 committed Jul 16, 2010 264 265 \end{lst}  Andreas Lauser committed Feb 13, 2012 266 267 First, the spatial parameters type tag is created on line \ref{tutorial-coupled:define-spatialparameters-typetag}. The type tag  Klaus Mosthaf committed Feb 23, 2012 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  Andreas Lauser committed Feb 13, 2012 270 271 parameters class itself (line \ref{tutorial-coupled:set-spatialparameters}) or the capillary  Klaus Mosthaf committed Feb 23, 2012 272 pressure/relative permeability relations\footnote{Taken together, the  Andreas Lauser committed Feb 13, 2012 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  Klaus Mosthaf committed Feb 23, 2012 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  Andreas Lauser committed Feb 13, 2012 286 requires, it provides a parameter class  Klaus Mosthaf committed Feb 23, 2012 287 288 \texttt{RegularizedBrooksCoreyParams} which has the type \texttt{Params} and which is defined in line  Andreas Lauser committed Feb 13, 2012 289 290 \ref{tutorial-coupled:matLawObjectType}. In this case, the spatial parameters only require a single set of parameters which means that it  Klaus Mosthaf committed Feb 23, 2012 291 only requires a single material parameter object as can be seen in  Andreas Lauser committed Feb 13, 2012 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  Klaus Mosthaf committed Feb 09, 2012 296 to be called with three arguments:  Andreas Lauser committed Feb 13, 2012 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  Klaus Mosthaf committed Feb 23, 2012 300  entity. Elements provide information about their geometry and  Andreas Lauser committed Feb 13, 2012 301  position and can be mapped to a global index.  Klaus Mosthaf committed Feb 23, 2012 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  Andreas Lauser committed Feb 13, 2012 307 308 309  the element. \end{description}  Klaus Mosthaf committed Feb 23, 2012 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  Andreas Lauser committed Feb 13, 2012 312 dimension of the grid.  Philipp Nuske committed Dec 10, 2010 313   Klaus Mosthaf committed Feb 09, 2012 314 The method \texttt{porosity()} defined in line  Bernd Flemisch committed Jul 16, 2010 315 \ref{tutorial-coupled:porosity} is called with the same arguments as  Andreas Lauser committed Feb 13, 2012 316 317 318 \texttt{intrinsicPermeability()} and returns a scalar value for porosity dependent on the position in the domain.  Klaus Mosthaf committed Feb 23, 2012 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  Andreas Lauser committed Feb 13, 2012 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  Benjamin Faigle committed Feb 04, 2011 334 \verb+dumux/material/fluidmatrixinteractions/2p+.  Bernd Flemisch committed Jul 16, 2010 335 336 337 338  \subsection{Exercises} \label{tutorial-coupled:exercises} The following exercises will give you the opportunity to learn how you  Andreas Lauser committed Feb 13, 2012 339 340 can change soil parameters, boundary conditions, run-time parameters and fluid properties in \Dumux.  Bernd Flemisch committed Jul 16, 2010 341 342  \subsubsection{Exercise 1}  Klaus Mosthaf committed Feb 09, 2012 343 344 \renewcommand{\labelenumi}{\alph{enumi})} For Exercise 1 you have to make only some small changes in the tutorial files.  Melanie Darcis committed Nov 11, 2010 345   Bernd Flemisch committed Jul 16, 2010 346 347 \begin{enumerate}  Philipp Nuske committed Feb 14, 2012 348 \item \textbf{Running the Program} \\  Christoph Grueninger committed Feb 24, 2012 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}.  Melanie Darcis committed Nov 11, 2010 350 351  \item \textbf{Changing the Model Domain and the Boundary Conditions} \\  Bernd Flemisch committed Jul 16, 2010 352  Change the size of the model domain so that you get a rectangle with  Philipp Nuske committed Feb 14, 2012 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 committed Jul 16, 2010 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  Klaus Mosthaf committed Feb 09, 2012 363  run the model as explained above.  Bernd Flemisch committed Jul 16, 2010 364   Klaus Mosthaf committed Feb 23, 2012 365  \item \textbf{Changing the Shape of the Discrete Elements} \\  Philipp Nuske committed Feb 14, 2012 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 committed Jul 16, 2010 368 \item \textbf{Changing Fluids} \\  Andreas Lauser committed Feb 13, 2012 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:  Benjamin Faigle committed Feb 04, 2011 370 \begin{enumerate}  Klaus Mosthaf committed Feb 23, 2012 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{}.  Andreas Lauser committed Feb 13, 2012 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.  Melanie Darcis committed Nov 11, 2010 373 \end{enumerate}  Viktor Szukitsch committed Sep 01, 2011 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  Andreas Lauser committed Feb 13, 2012 375 \texttt{/dumux/material/components} and read chapter~\ref{sec:fluidframework}.  Bernd Flemisch committed Jul 16, 2010 376   Andreas Lauser committed Feb 13, 2012 377 \item \textbf{Use a Full-Fledged Fluid System} \\  Andreas Lauser committed Feb 13, 2012 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.\\  Klaus Mosthaf committed Feb 23, 2012 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}. 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.\\  Andreas Lauser committed Feb 13, 2012 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 committed Jul 16, 2010 381 382  \item \textbf{Changing Constitutive Relations} \\  Andreas Lauser committed Feb 13, 2012 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  Philipp Nuske committed Dec 10, 2010 384  regularized Brooks-Corey law for the  Melanie Darcis committed Nov 11, 2010 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 committed Jul 16, 2010 387  You can find the material laws in the folder  Melanie Darcis committed Oct 13, 2010 388  \verb+dumux/material/fluidmatrixinteractions+. The necessary parameters  Melanie Darcis committed Nov 11, 2010 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 committed Jul 16, 2010 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  Melanie Darcis committed Nov 11, 2010 396  so that water is again flowing from the left to the right of the  397 \begin{figure}[ht]  Andreas Lauser committed Feb 13, 2012 398 \psfrag{K1 =}{$\mathbf{K} = 10^{-8}\;\text{m}^2$}  Bernd Flemisch committed Jul 16, 2010 399 \psfrag{phi1 =}{$\phi = 0.15$}  Andreas Lauser committed Feb 13, 2012 400 \psfrag{K2 =}{\textcolor{white}{$\mathbf{K} = 10^{-9}\;\text{m}^2$}}  Bernd Flemisch committed Jul 16, 2010 401 \psfrag{phi2 =}{\textcolor{white}{$\phi = 0.3$}}  Andreas Lauser committed Feb 13, 2012 402 403 \psfrag{600 m}{$600 \;\text{m}$} \psfrag{300 m}{$300 \;\text{m}$}  Bernd Flemisch committed Jul 16, 2010 404 405 \centering \includegraphics[width=0.5\linewidth,keepaspectratio]{EPS/exercise1_c.eps}  Andreas Lauser committed Feb 13, 2012 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 committed Jul 16, 2010 407 \end{figure}  Benjamin Faigle committed Feb 04, 2011 408 domain. You can use the fluids of exercise 1c).\\  Andreas Lauser committed Feb 13, 2012 409 \textbf{Hint:} The current position of the control volume can be obtained using \texttt{element\allowbreak.geometry()\allowbreak.corner(scvIdx)}.\\  Andreas Lauser committed Feb 13, 2012 410 When does the front cross the material border? In paraview, the  Andreas Lauser committed Feb 13, 2012 411 animation view (\textit{View} $\rightarrow$ \textit{Animation  Andreas Lauser committed Feb 13, 2012 412 413  View}) is a convenient way to get a rough feeling of the time-step sizes.  Bernd Flemisch committed Jul 16, 2010 414 415 416 \end{enumerate} \subsubsection{Exercise 2}  Klaus Mosthaf committed Feb 25, 2011 417 For this exercise you should create a new problem file analogous to  Andreas Lauser committed Feb 13, 2012 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  Andreas Lauser committed Feb 13, 2012 421 problem file needs to  Klaus Mosthaf committed Feb 23, 2012 422 423 be included in the file \texttt{tutorial\_coupled.cc}.  Andreas Lauser committed Feb 13, 2012 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}  Klaus Mosthaf committed Feb 23, 2012 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  Andreas Lauser committed Feb 13, 2012 433 use a new type tag for the problem as well as a new problem class  Klaus Mosthaf committed Feb 23, 2012 434 e.g. \mbox{\texttt{Ex2TutorialProblemCoupled}}. Make sure to assign your  Andreas Lauser committed Feb 13, 2012 435 newly defined spatial parameter class to the  436 \texttt{SpatialParams} property for the new  Klaus Mosthaf committed Feb 23, 2012 437 438 type tag.  Andreas Lauser committed Feb 13, 2012 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 committed Jul 16, 2010 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]  Andreas Lauser committed Feb 13, 2012 455 \psfrag{K1}{K $= 10^{-7}\;\text{m}^2$}  Bernd Flemisch committed Jul 16, 2010 456 \psfrag{phi1}{$\phi = 0.2$}  Klaus Mosthaf committed Feb 23, 2012 457 \psfrag{Lin}{\textsc{Brooks}-\textsc{Corey} Law}  Andreas Lauser committed Feb 13, 2012 458 459 \psfrag{Lin2}{$\lambda = 1.8$, $p_e = 1000\;\text{Pa}$} \psfrag{K2}{K $= 10^{-9}\;\text{m}^2$}  Bernd Flemisch committed Jul 16, 2010 460 \psfrag{phi2}{$\phi = 0.15$}  Andreas Lauser committed Feb 13, 2012 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 committed Jul 16, 2010 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]  Andreas Lauser committed Feb 13, 2012 475 \psfrag{pw}{$p_w = 5 \times 10^5\;\text{Pa}$}  Bernd Flemisch committed Jul 16, 2010 476 \psfrag{S}{$S_n = 1.0$}  Andreas Lauser committed Feb 13, 2012 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 committed Jul 16, 2010 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}  Melanie Darcis committed Nov 11, 2010 485 \begin{itemize}  Andreas Lauser committed Feb 13, 2012 486  \item Increase the simulation time to e.g. $4\times 10^7 \;\text{s}$. Investigate the saturation: Is the value range reasonable?  Melanie Darcis committed Nov 11, 2010 487 488 489  \item What happens if you increase the resolution of the grid? \end{itemize}  Philipp Nuske committed Feb 14, 2012 490 \subsubsection{Exercise 3: Parameter File Input}  Andreas Lauser committed Feb 13, 2012 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 committed Feb 24, 2012 494 via \textit{parameter input files}.  Benjamin Faigle committed Feb 12, 2012 495   Andreas Lauser committed Feb 13, 2012 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  Philipp Nuske committed Feb 14, 2012 503 happens if they are modified. For more information about the input file please refer to section \ref{sec:inputFiles}.  Benjamin Faigle committed Feb 12, 2012 504   Philipp Nuske committed Feb 14, 2012 505 \subsubsection{Exercise 4: Create a New Component}  Bernd Flemisch committed Jul 16, 2010 506   Andreas Lauser committed Feb 13, 2012 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  Klaus Mosthaf committed Feb 23, 2012 509 existing components in the directory \verb+/dumux/material/components+). \\  Bernd Flemisch committed Jul 16, 2010 510 Use benzene as a new fluid and run the model of Exercise 2 with water  Andreas Lauser committed Feb 13, 2012 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 committed Jul 16, 2010 513   514 \clearpage \newpage  515 516 517 518 %%% Local Variables: %%% mode: latex %%% TeX-master: "dumux-handbook" %%% End: