opm-simulators/doc/handbook/tutorial-decoupled.tex

\clearpage
\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}

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 as well as capillarity effects are neglected.

\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]$}
\psfrag{q_n = -3 x 10^-4 [kg/m^2s]}{$q_n = -3 \times 10^-4$ $\left[\frac{\textnormal{kg}}{\textnormal{m}^2 \textnormal{s}}\right]$}
\centering
\includegraphics[width=0.9\linewidth,keepaspectratio]{EPS/tutorial-problemconfiguration}
\caption{Geometry of the tutorial problem with initial and boundary conditions.}\label{tutorial-decoupled:problemfigure}
\end{figure}

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{/dune-mux/test/tutorial}.

\begin{lst}[File dune-mux/test/tutorial/tutorial\_decoupled.cc]\label{tutorial-deoucpled:mainfile} \mbox{}
\lstinputlisting[basicstyle=\ttfamily\scriptsize,numbers=left, 
numberstyle=\tiny, numbersep=5pt]{../../test/tutorial/tutorial_decoupled.cc}
\end{lst}

First, in line \ref{tutorial-decoupled:include-begin} to \ref{tutorial-decoupled:include-end} the Dune and the \Dumux files respectively which contain the functions and classes that are needed in the main loop are included into the main file.

In line \ref{tutorial-decoupled:grid-begin} to \ref{tutorial-decoupled:grid-end} 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 of cells in $x$ and $y$ direction (\texttt{N}), where the dimensions are previously defined in line \ref{tutorial-decoupled:dim}. The grid of type \texttt{Dune::SGrid} is then generated in line \ref{tutorial-decoupled:grid-end}. For more information about the dune grid interface, the different grid types that are supported and the generation of different grids it is referred to the \textit{Dune Grid Interface HOWTO} \cite{DUNE-HP}.

The second point mentioned at the beginning of this section was the definition of material properties and constitutive relationships. The fluid properties of the two fluid phases considered here are defined in lines \ref{tutorial-decoupled:water} and \ref{tutorial-decoupled:oil}. The fluid classes including different kinds of fluids (here: \texttt{Dune::Water} and \texttt{Dune::Oil}) can be found in the file \texttt{phaseproperties2p.hh} in the directory \texttt{/dune-mux/dumux}-\texttt{/material/phaseproperties}. The properties of the solid matrix are defined in a special soil class. The \texttt{soil} object is generated in line \ref{tutorial-decoupled:soil}. As can be seen, the class type is \texttt{Dune::TutorialSoil}, which is defined in the file \texttt{tutorial\_soilproperties\_decoupled.hh} in the folder \texttt{/test/tutorial}. A description of this file and the definition of a soil class including the soil parameters can be found in section \ref{tutorial-decoupled:description-soil-class}. Finally, in line \ref{tutorial-decoupled:twophaserelations} the information included in the fluid and soil objects is used to generate an object of type \texttt{Dune::TwoPhaseRelations}, which includes the constitutive relationships (capillary pressure-saturation relation, relative permeability-saturation relation, etc.). The file \texttt{twophaserelations.hh} can be found in the directory \texttt{/dune-mux/dumux/material}.

The definition of boundary and initial conditions as well as source of sink terms is done by definition of a so-called \textit{problem} class. In case of this tutorial the problem class is defined in the file \texttt{tutorialproblem\_decoupled.hh} in the \texttt{/test/tutorial} folder. In the main file the problem object of type \texttt{Dune::TutorialProblemDecoupled} is then generated in line \ref{tutorial-decoupled:problem}. A further explanation of the definition of boundary and initial conditions, source and sink terms and the structure of the problem class can be found in section \ref{tutorial-decoupled:description-bc-ic}. Besides the definition of the boundary and initial conditions the problem class is also a kind of interface containing all the objects generated before (geometry, fluids, soil, constitutive relationships, etc.). Thus, as can be seen in line \ref{tutorial-decoupled:problem} all this objects are given as arguments when calling the constructor of the problem class.

Following the steps listed at the beginning of this section, finally, a model has to be chosen. In case of this tutorial a decoupled isothermal two phase model is the choice. As explained before a decoupled model consists of a pressure equation which is decoupled or only weakly coupled to a saturation equation, concentration equations, energy balance equations, etc. In case of isothermal two phase flow one pressure equation and one saturation equation have to be solved. 

The discretisation of the pressure equation is included in the object which is generated in line \ref{tutorial-decoupled:diffusion} of the main file. It is called \texttt{diffusion} and it is of type \texttt{Dune::FVTotalVelocity2P}. The definition of this class can be found in \texttt{/dune-mux/dumux/diffusion/fv} in the file \texttt{fvtotalvelocity2p.hh}. The \texttt{FV} in the class name indicates that a \textit{finite volume} discretisation is used. For the pressure equation also other discretisation methods like \textit{finite elements} or \textit{mimetic finite differences} are available.

Similarly, an object containing the discretisation of the saturation equation is generated in line \ref{tutorial-decoupled:transport} of the main file. As can be seen, this object is called \texttt{transport}. The class  \texttt{Dune::FVSaturationWetting2P} can be found in the directory \texttt{/dune-mux/dumux/transport/fv}. Again, the \texttt{FV} in the class name indicates a \textit{finite volume} discretisation.

For the decoupled model the so-called \textit{IMPES} (\textit{IM}plicit \textit{P}ressure \textit{E}xplicit \textit{S}aturation) algorithm is used. This means that the pressure equation is first solved using an implicit method. The resulting velocities are then used to solve the saturation equation explicitly. It is clear that depending on the non-linearity this is an iterative process within every timestep. The object \texttt{impes} of type \texttt{Dune::IMPES} which is generated in line \ref{tutorial-decoupled:impes} of the main file combines the objects \texttt{diffusion} and \texttt{transport} and proceeds the iterations of the IMPES algorithm within a timestep. The definition of the class \texttt{Dune::IMPES} can be found in the file \texttt{impes.hh} in the directory \texttt{/dune-mux/dumux/fractionalflow/impes}.

Finally, an object called \texttt{timeloop} of type \texttt{Dune::TimeLoop} is generated in line \ref{tutorial-decoupled:timeloop} of the tutorial main file. The class \texttt{Dune::TimeLoop} is defined in the file \texttt{timeloop.hh} in the folder \texttt{/dune-mux/dumux/timedisc}. The object \texttt{timeloop} includes the type of timestep that is used (implicit, explicit, etc.) and contains the function \texttt{execute} which is called in line \ref{tutorial-decoupled:execute} of the main file. This function finally starts the computation and runs the (time)loop over all timesteps.

\subsection{The definition of the fluid properties}\label{tutorial-decoupled:description-fluid-class}

In \Dumux different fluids are already implemented. The definitions can be found in the file \texttt{phaseproperties2p.hh} in the directory \texttt{/dune-mux/dumux/material/phaseproperties}. As can be observed in this file, for each fluid a class named like the fluid is defined. These classes are derived from the fluid base class \texttt{Fluid} which is defined in the file \texttt{property\_baseclasses.hh} in the directory \texttt{/dune-mux/dumux/material} and include several functions returning different fluid properties.

If a new fluid type is to be used or an existing type is to be used with different fluid parameters new fluid types can be defined derived from the base class \texttt{Fluid} comparable to the fluids already defined in the file \texttt{phaseproperties2p.hh}.

It is important to mention, that existing fluid classes should not be changed. New fluid classes should only be added to the file \texttt{phaseproperties2p.hh} if they are also to be added to the repository! If you are not sure if your fluid class can be useful for the other \Dumux users just create a new file in your problem directory similar to the file \texttt{phaseproperties2p.hh} and define your fluid classes there. 

\subsection{The definition of the soil parameters}\label{tutorial-decoupled:description-soil-class}

Soil properties which can be defined in \Dumux are the \textit{intrinsic permeability}, the \textit{porosity} and the \textit{heat capacity} as well as the \textit{heat conductivity} of the solid matrix. Further the \textit{residual saturations} of the fluids, and the \textit{capillary pressures-saturation function} as well as the \textit{relative permeability-saturation functions} are depending on the soil.

The base class \texttt{Dune::Matrix2p} for the definition of the soil parameters can be found in the file \texttt{property\_baseclasses.hh} in the directory \texttt{/dune-mux/dumux/material}. Derived from this base class, there exist two standard soil type classes named \texttt{HomogeneousSoil} and \texttt{HeterogeneousSoil}. Both can be found in the file \texttt{matrixproperties.hh} in the \texttt{/material} folder. If one wants to use a soil that differs from this standard soil types, new soil classes can be derived either from the base class (\texttt{Dune::Matrix2p}) or from the two standard soil classes (\texttt{Dune::HomogeneousSoil} and \texttt{Dune::HeterogeneousSoil}).

For this tutorial problem a new soil class named \texttt{TutorialSoil} is derived from \texttt{Dune::HomogeneousSoil} (listing \ref{tutorial-deoucpled:soilpropertiesfile}, line \ref{tutorial-decoupled:tutorialsoil}), which can be found in the file \texttt{tutorial\_soilproperties\_decoupled.hh} in the directory \texttt{/test/tutorial}.

Listing \ref{tutorial-deoucpled:soilpropertiesfile} shows the file \texttt{tutorial\_soilproperties\_decoupled.hh}.

\begin{lst}[File dune-mux/test/tutorial/tutorial\_soilproperties\_decoupled.hh]\label{tutorial-deoucpled:soilpropertiesfile} \mbox{}
\lstinputlisting[basicstyle=\ttfamily\scriptsize,numbers=left, 
numberstyle=\tiny, numbersep=5pt]{../../test/tutorial/tutorial_soilproperties_decoupled.hh}
\end{lst}

In line \ref{tutorial-decoupled:permeability} the function returning the intrinsic permeability can be found. As can be seen, the function has to be called with three different arguments. The first one (\texttt{x}) is a vector including the global coordinates of the current entity (can be an element, vertex, etc.), the second one (\texttt{e}) is the entity itself and the third one is a vector including the local coordinates of the current entity. The intrinsic permeability is a tensor and thus returned in form of a $n \times n$-matrix where $n$ is the dimension of the problem.

The function \texttt{porosity()} defined in line \ref{tutorial-decoupled:porosity} is called with the same arguments as the permeability function described before and returns the porosity dependent on the position in the domain.

The residual saturation functions \texttt{Sr\_w()} (line \ref{tutorial-decoupled:srw}) and \texttt{Sr\_n()} (line \ref{tutorial-decoupled:srn}) additionally have the temperature as function argument, which is set to a default value if a isothermal model is used.

Finally, the functions defining the type of the capillary pressure function and the relative permeability functions have to be considered. In line \ref{tutorial-decoupled:flags} the function \texttt{relPermFlag()} is defined. This function returns a flag indicating the type of function which is used depending on the position. This could be a linear function, a \textit{Brooks-Corey} function, a \textit{van Genuchten} function, etc. The flags that can be chosen as return parameter are defined in the base soil class \texttt{Matrix2p} in the file \texttt{property\_baseclasses.hh}. The parameters used in the chosen function type can be defined in the function \texttt{paramRelPerm} (line \ref{tutorial-decoupled:parameters}). As can be seen in listing \ref{tutorial-deoucpled:soilpropertiesfile}, e.g. linear capillary pressure and relative permeability functions require a vector of two arguments, one defining the minimum and one defining the maximum capillary pressure. The parameters can again be defined depending on the position in the domain an on temperature.

\subsection{The definition of boundary and initial conditions and source or sink terms}\label{tutorial-decoupled:description-bc-ic}

Boundary and initial conditions are defined in a so-called problem class. The problem class of this tutorial has the name \texttt{TutorialProblemDecoupled} and is defined in the file \texttt{tutorialproblem\_decoupled.hh} which can be found in the directory \texttt{/test/tutorial}. Listing \ref{tutorial-deoucpled:problemfile} shows the class \texttt{TutorialProblemDecoupled}. As can be seen it is derived from the problem base class \texttt{FractionalFlowProblem} (line \ref{tutorial-decoupled:tutorialproblem}) which is defined in the file \texttt{fractionalflowproblem.hh} in the directory \texttt{/dune-mux/dumux/fractionalflow}. 

\begin{lst}[File dune-mux/test/tutorial/tutorialproblem\_decoupled.hh]\label{tutorial-deoucpled:problemfile} \mbox{}
\lstinputlisting[basicstyle=\ttfamily\scriptsize,numbers=left, 
numberstyle=\tiny, numbersep=5pt]{../../test/tutorial/tutorialproblem_decoupled.hh}
\end{lst}

Listing \ref{tutorial-decoupled:tutorialproblem}) includes five types of functions. The type of each function can be identified by the first letter and the first part respectively of its name. Function names of functions returning
\begin{itemize}
 \item a source or sink term start with a \textbf{q},
 \item a boundary condition type start with \textbf{bctype},
 \item a \textit{Dirichlet} boundary condition start with a \textbf{g},
 \item a \textit{Neumann} boundary condition start with a \textbf{J} and
 \item a initial condition start with \textbf{init}. 
\end{itemize}
The second part of the function name indicates if the return value is used for the pressure equation (\textbf{Press}) or for the saturation equation (\textbf{Sat}).

All different function types have to be called with three different arguments. The first one (\texttt{x}) is a vector including the global coordinates of the current entity (can be an element, vertex, etc.), the second one (\texttt{e}) is the entity itself and the third one is a vector including the local coordinates of the current entity. Thus, the return of the functions, which can be a boundary value, an initial value, a source/sink, etc., can be defined depending on the position in the domain.

The first function defined in the problem class \texttt{TutorialProblemDecoupled} is the function \texttt{qPress} (line \ref{tutorial-decoupled:qpress}). It returns a source or a sink term for the pressure equation. 

In lines \ref{tutorial-decoupled:bctypepress} and \ref{tutorial-decoupled:bctypesat} the functions returning the boundary condition type are defined. Flags of type \texttt{Dune::BoundaryConditions::Flags} have to be used as return value of these functions. The flags that can be chosen are defined in the file \texttt{boundaryconditions.hh} in the directory \texttt{/dune-disc/disc/operators}.

In lines \ref{tutorial-decoupled:gpress} and \ref{tutorial-decoupled:gsat} the functions returning the \textit{Dirichlet} boundary conditions and in line \ref{tutorial-decoupled:jpress} a function returning the \textit{Neumann} boundary conditions are defined.

Finally, the function \texttt{initSat} is defined in line \ref{tutorial-decoupled:initsat}. This function returns the initial saturation distribution.

\subsection{Exercise}
\label{tutorial-deoucpled:exercises}
The following exercises will give you the opportunity to learn how you can change soil parameters, boundary conditions and fluid properties in \Dumux. For each exercise you can find the output file of the last timestep in the directory \texttt{/dune-mux/dumux/tutorial/results/decoupled}.

\subsubsection{Exercise 1}
\renewcommand{\labelenumi}{\alph{enumi})}
For Exercise 1 you only have to make some small changes in the tutorial files.
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}. For the visualisation with paraview please refer to \ref{quick-start-guide}.
\begin{enumerate}
\item \textbf{Changing the Model Domain and the Boundary Conditions} \\
Change the size of the model domain so that you get a rectangle
with edge lengths of x = 400 m \\  and y = 500 m and with discretisation lengths of  $\Delta \text{x} = 20$ m and $\Delta \text{y} = 20$ m. \\
Change the boundary conditions in the file \texttt{tutorialproblem\_decoupled.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\_decoupled} and run the model.

\item \textbf{Changing Fluids} \\
Now you can change the fluids. Use DNAPL instead of Oil and Brine instead of Water. To do that you have to change the file \texttt{tutorial\_decoupled.cc}. If you want to take a closer look how the fluid classes are defined and which fluids are already available please open the file \texttt{phaseproperties2p.hh} in the directory
\texttt{/dune-mux/dumux/material/phaseproperties}.

\item \textbf{Changing Constitutive Relationships} \\
Use a Brooks-Corey law with $\lambda$ = 2 and entry pressure $p_b = 0.0$ instead of a linear law for the relative-permeability/saturation relationship. To do that you have to change the file \texttt{tutorial\_soilproperties\_decoupled.hh}. You can find the flag that you have to set for the Brooks-Corey law in the file \texttt{property\_baseclasses.hh} in the directory \texttt{/dune-mux/dumux/material}.  
The available relative permeability and capillary pressure functions are defined in the file \texttt{/dune-mux/dumux/material/relperm\_pc\_law}.\\
(Hint: The CFL-security-factor should be set to a value smaller than 1 if nonlinear relationships are used (Suggestion: 0.95). It can be found in the application file \texttt{tutorialproblem\_decoupled.hh}.)
 
\item \textbf{Heterogeneities}  \\
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 still flowing from the bottom to the top of the domain. You can use the fluids of exercise 1b) and the constitutive relationship of exercise 1c). 

\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}

\end{enumerate}

\subsubsection{Exercise 2}
For this exercise you should create a new proplem file according to the file \texttt{tutorialproblem\_decoupled.hh} and a new soil property file according to the file \texttt{tutorial\_soilproperties\_decoupled.hh}. These files need to be included in the file \texttt{tutorial\_decoupled.cc}.\\
The new soil file should contain the definition of a new soil class e.g. SoilEx2. Make sure that you also adjust the preprocessor commands (e.g. change TUTORIAL\_SOILPROPERTIES to TUTORIAL\_SOILEX2). 
The new problem file should contain the definition of a new problem class e.g. ProblemEx2. Here you also need to adjust the preprocessor commands. 
Replace the classes \texttt{TutorialSoil} and \texttt{TutorialProblemDecoupled} by the the new classes you just created. \\
Now, set up a model that describes the processes given in the Figures \ref{tutorial-deoucpled:ex2_Domain} and \ref{tutorial-deoucpled:ex2_BC}. Initially the domain is fully saturated with water and the pressure is $p_w = 5 \times 10^5$ 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 $4\times 10^7$ s.\\
(Hint: Chose the CFL-security-factor as 0.95 - see Exercise 1 c).)

\begin{figure}[h]
\psfrag{K1}{K $= 10^{-7}\text{ m}^2$}
\psfrag{phi1}{$\phi = 0.2$}
\psfrag{Lin}{Linear Law}
\psfrag{K2}{K $= 10^{-9}\text{ m}^2$}
\psfrag{phi2}{$\phi = 0.15$}
\psfrag{BC1}{Brooks Corey Law} 
\psfrag{BC2}{$\lambda = 1.8$, $p_b = 0.0$}
\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}
\caption{Set-up of the model domain and the soil parameters}\label{tutorial-deoucpled:ex2_Domain}
\end{figure}

\begin{figure}[h]
\psfrag{pw}{$p_w = 5 \times 10^5$ \text{Pa}}
\psfrag{S}{$S_n = 1.0$}
\psfrag{qw}{$q_w = 2 \times 10^{-7}$ [kg/$\text{m}^2$s]}
\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}
\caption{Boundary Conditions}\label{tutorial-deoucpled:ex2_BC}
\end{figure}

\subsubsection{Exercise 3}
Create a file called \texttt{new\_fluid.hh} and implement a new fluid class. This new fluid class should be derived from the base class Fluid which can be found in  \texttt{/dune-mux/dumux/material/property\_baseclasses.hh}. \\
(You can look at existing fluid classes in the file \texttt{/dune-mux/dumux/material/phaseproperties/phaseproperties2p.hh}.)
Use benzene as a new fluid and run the model of Exercise 2 with water and benzene. The properties of benzene are given in the following: \\
density: 889.51 kg/$\text{m}^3$, viscosity: 0.00112 Pa s.