diff --git a/doc/handbook/3_tutorialimplicit.tex b/doc/handbook/3_tutorialimplicit.tex index f908bceec97881751f9a8f23ec3efc04ada86a40..387871652d2d8aa61aad34b03171635483f9fb7c 100644 --- a/doc/handbook/3_tutorialimplicit.tex +++ b/doc/handbook/3_tutorialimplicit.tex @@ -65,12 +65,12 @@ The solved equations are the mass balances of water and oil: \subsection{The Main File} Listing \ref{tutorial-implicit:mainfile} shows the main application file -\texttt{tutorial/tutorial\_implicit.cc} for the implicit two-phase +\texttt{tutorial/tutorial\_implicit/tutorial\_implicit.cc} for the implicit two-phase model. This file has to be compiled and executed in order to solve the problem described above. -\begin{lst}[File tutorial/tutorial\_implicit.cc]\label{tutorial-implicit:mainfile} \mbox{} - \lstinputlisting[style=DumuxCode, numbersep=5pt, firstline=24, firstnumber=24]{../../tutorial/tutorial_implicit.cc} +\begin{lst}[File tutorial/tutorial\_implicit/tutorial\_implicit.cc]\label{tutorial-implicit:mainfile} \mbox{} + \lstinputlisting[style=DumuxCode, numbersep=5pt, firstline=24, firstnumber=24]{../../tutorial/tutorial_implicit/tutorial_implicit.cc} \end{lst} From line \ref{tutorial-implicit:include-begin} to line @@ -99,8 +99,8 @@ values provided in the command line have precedence. Listing~\ref{tutorial-implicit:parameter-file} shows the default parameter file for the tutorial problem. -\begin{lst}[File tutorial/tutorial\_implicit.input]\label{tutorial-implicit:parameter-file} \mbox{} -\lstinputlisting[style=DumuxParameterFile]{../../tutorial/tutorial_implicit.input} +\begin{lst}[File tutorial/tutorial\_implicit/tutorial\_implicit.input]\label{tutorial-implicit:parameter-file} \mbox{} +\lstinputlisting[style=DumuxParameterFile]{../../tutorial/tutorial_implicit/tutorial_implicit.input} \end{lst} To provide an error message, the usage message which is displayed to @@ -116,8 +116,8 @@ When solving a problem using \Dumux, the most important file is the so-called \textit{problem file} as shown in listing~\ref{tutorial-implicit:problemfile}. -\begin{lst}[File tutorial/tutorialproblem\_implicit.hh]\label{tutorial-implicit:problemfile} \mbox{} -\lstinputlisting[style=DumuxCode, numbersep=5pt, firstline=24, firstnumber=24]{../../tutorial/tutorialproblem_implicit.hh} +\begin{lst}[File tutorial/tutorial\_implicit/tutorialproblem\_implicit.hh]\label{tutorial-implicit:problemfile} \mbox{} +\lstinputlisting[style=DumuxCode, numbersep=5pt, firstline=24, firstnumber=24]{../../tutorial/tutorial_implicit/tutorialproblem_implicit.hh} \end{lst} First, a new type tag is created for the problem in line @@ -247,7 +247,7 @@ 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 +\Dumux fluid framework can be found at \url{http://www.dumux.org/doxygen-stable/html-\DumuxVersion/modules.php} %in chapter~\ref{sec:fluidframework}. @@ -270,8 +270,8 @@ should be derived from the base class \ref{tutorial-implicit:spatialparametersfile} shows the file \\ \verb+tutorialspatialparams_implicit.hh+: -\begin{lst}[File tutorial/tutorialspatialparams\_implicit.hh]\label{tutorial-implicit:spatialparametersfile} \mbox{} -\lstinputlisting[style=DumuxCode, numbersep=5pt, firstline=25, firstnumber=25]{../../tutorial/tutorialspatialparams_implicit.hh} +\begin{lst}[File tutorial/tutorial\_implicit/tutorialspatialparams\_implicit.hh]\label{tutorial-implicit:spatialparametersfile} \mbox{} +\lstinputlisting[style=DumuxCode, numbersep=5pt, firstline=25, firstnumber=25]{../../tutorial/tutorial_implicit/tutorialspatialparams_implicit.hh} \end{lst} First, the spatial parameters type tag is created on line @@ -411,7 +411,7 @@ To do that, you have to select different components via the property system in t \end{enumerate} 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 -\texttt{/dumux/material/components} and read +\texttt{/dumux/material/components} and read the doxygen documentation \url{http://www.dumux.org/doxygen-stable/html-\DumuxVersion/modules.php}. %chapter~\ref{sec:fluidframework}. @@ -597,7 +597,7 @@ As you have experienced, compilation takes quite some time. Therefore, via \textit{parameter input files}. In the code, parameters can be read via the macro -\texttt{GET\_RUNTIME\_PARAM(TypeTag, Scalar, +\texttt{GET\_RUNTIME\_PARAM(TypeTag, Scalar, MyWonderful\allowbreak Group.MyWonderfulParameter);}. In this exercise we will explore the possibilities of the parameter file. For this we take a look at the file \texttt{ex3\_tutorial\_implicit.input} in the \texttt{solutions\_implicit} folder. Besides the parameters which you already used in the parameter file above, diff --git a/doc/handbook/3_tutorialsequential.tex b/doc/handbook/3_tutorialsequential.tex index e2d8754671ce581cc02285e7db6c1d3a750e62ea..228ddb3bc9f6fa2b4b997d04b5e6d7adc5e152d1 100644 --- a/doc/handbook/3_tutorialsequential.tex +++ b/doc/handbook/3_tutorialsequential.tex @@ -57,8 +57,8 @@ executed, has to be set up, if the problem described above is to be solved using a sequential model. This main file can be found in the directory \texttt{/tutorial} of the stable part of \Dumux. -\begin{lst}[File tutorial/tutorial\_sequential.cc]\label{tutorial-sequential:mainfile} \mbox{} -\lstinputlisting[style=DumuxCode, numbersep=5pt, firstline=24, firstnumber=24]{../../tutorial/tutorial_sequential.cc} +\begin{lst}[File tutorial/tutorial\_sequential/tutorial\_sequential.cc]\label{tutorial-sequential:mainfile} \mbox{} +\lstinputlisting[style=DumuxCode, numbersep=5pt, firstline=24, firstnumber=24]{../../tutorial/tutorial_sequential/tutorial_sequential.cc} \end{lst} First, from line \ref{tutorial-sequential:include-begin} to line @@ -85,8 +85,8 @@ values provided in the command line have precedence. Listing~\ref{tutorial-sequential:parameter-file} shows the default parameter file for the tutorial problem. -\begin{lst}[File tutorial/tutorial\_sequential.input]\label{tutorial-sequential:parameter-file} \mbox{} -\lstinputlisting[style=DumuxParameterFile]{../../tutorial/tutorial_sequential.input} +\begin{lst}[File tutorial/tutorial\_sequential/tutorial\_sequential.input]\label{tutorial-sequential:parameter-file} \mbox{} +\lstinputlisting[style=DumuxParameterFile]{../../tutorial/tutorial_sequential/tutorial_sequential.input} \end{lst} To provide an error message, the usage message which is displayed to @@ -104,8 +104,8 @@ so-called \textit{problem file} as shown in listing \ref{tutorial-sequential:problemfile} of \texttt{tutorialproblem\_sequential.hh}. -\begin{lst}[File tutorial/tutorialproblem\_sequential.hh]\label{tutorial-sequential:problemfile} \mbox{} -\lstinputlisting[style=DumuxCode, numbersep=5pt, firstline=24, firstnumber=24]{../../tutorial/tutorialproblem_sequential.hh} +\begin{lst}[File tutorial/tutorial\_sequential/tutorialproblem\_sequential.hh]\label{tutorial-sequential:problemfile} \mbox{} +\lstinputlisting[style=DumuxCode, numbersep=5pt, firstline=24, firstnumber=24]{../../tutorial/tutorial_sequential/tutorialproblem_sequential.hh} \end{lst} First, both \Dune grid handlers and the sequential model of \Dumux @@ -189,8 +189,8 @@ on the functions, consult the documentation in the code. Listing \ref{tutorial-sequential:spatialparamsfile} shows the file \verb+tutorialspatialparams_sequential.hh+: -\begin{lst}[File tutorial/tutorialspatialparams\_sequential.hh]\label{tutorial-sequential:spatialparamsfile} \mbox{} -\lstinputlisting[style=DumuxCode, numbersep=5pt, firstline=24, firstnumber=24]{../../tutorial/tutorialspatialparams_sequential.hh} +\begin{lst}[File tutorial/tutorial\_sequential/tutorialspatialparams\_sequential.hh]\label{tutorial-sequential:spatialparamsfile} \mbox{} +\lstinputlisting[style=DumuxCode, numbersep=5pt, firstline=24, firstnumber=24]{../../tutorial/tutorial_sequential/tutorialspatialparams_sequential.hh} \end{lst} As this file only slightly differs from the implicit version, consult chapter \ref{tutorial-implicit:description-spatialParameters} for explanations. @@ -211,18 +211,18 @@ For Exercise 1 you only have to make some small changes in the tutorial files. \item \textbf{Altering output} To get an impression what the results should look like you can first run the -original version of the sequential tutorial model by typing \texttt{./tutorial\_sequential}. +original version of the sequential tutorial model by typing \texttt{./tutorial/tutorial\_sequential/tutorial\_sequential}. The runtime parameters which are set can be found in the input file (listing~\ref{tutorial-sequential:parameter-file}). -If the input file has the same name than the main file (e.g. \texttt{tutorial\_sequential.cc} -and \texttt{tutorial\_sequential.input}), it is automatically chosen. If the name differs -the program has to be started typing \texttt{./tutorial\_sequential -parameterFile .input}. -For more options you can also type \texttt{./tutorial\_sequential -h}. For the +If the input file has the same name than the main file (e.g. \texttt{tutorial/tutorial\_sequential/tutorial\_sequential.cc} +and \texttt{tutorial/tutorial\_sequential/tutorial\_sequential.input}), it is automatically chosen. If the name differs +the program has to be started typing \texttt{./tutorial/tutorial\_sequential/tutorial\_sequential -parameterFile .input}. +For more options you can also type \texttt{./tutorial/tutorial\_sequential/tutorial\_sequential -h}. For the visualisation with paraview please refer to \ref{quick-start-guide}.\\ As you can see, the simulation creates many output files. To reduce these in order to perform longer simulations, change the method responsible for output (line \ref{tutorial-sequential:outputinterval} in the file \texttt{tutorialproblem\_\allowbreak sequential}) as to write an output only every 20 time-steps. Compile the main file by typing -\texttt{make tutorial\_sequential} and run the model. Now, run the simulation for 5e5 seconds. +\texttt{make tutorial/tutorial\_sequential/tutorial\_sequential} and run the model. Now, run the simulation for 5e5 seconds. \item \textbf{Changing the Model Domain and the Boundary Conditions} \\ Change the size of the model domain so that you get a rectangle @@ -527,4 +527,3 @@ so both models can also work with functions based ond \mbox{\texttt{...AtPos(Glo no matter if we model implicit or sequential. Try to formulate a spatial parameters file that works with both problems, the implicit and the sequential. Therein, only use functions at the position. - diff --git a/doc/handbook/CMakeLists.txt b/doc/handbook/CMakeLists.txt index f2f89dcb575b7b07f9ee84414275c978b9f7da04..e19249dee114cdd24f49d9a733270c02ed434494 100644 --- a/doc/handbook/CMakeLists.txt +++ b/doc/handbook/CMakeLists.txt @@ -22,14 +22,14 @@ set(TEX_INPUTS 5_models.tex 5_propertysystem.tex 5_spatialdiscretizations.tex - ../../tutorial/tutorial_implicit.cc - ../../tutorial/tutorial_implicit.input - ../../tutorial/tutorialproblem_implicit.hh - ../../tutorial/tutorialspatialparams_implicit.hh - ../../tutorial/tutorial_sequential.cc - ../../tutorial/tutorial_sequential.input - ../../tutorial/tutorialproblem_sequential.hh - ../../tutorial/tutorialspatialparams_sequential.hh) + ../../tutorial/tutorial_implicit/tutorial_implicit.cc + ../../tutorial/tutorial_implicit/tutorial_implicit.input + ../../tutorial/tutorial_implicit/tutorialproblem_implicit.hh + ../../tutorial/tutorial_implicit/tutorialspatialparams_implicit.hh + ../../tutorial/tutorial_sequential/tutorial_sequential.cc + ../../tutorial/tutorial_sequential/tutorial_sequential.input + ../../tutorial/tutorial_sequential/tutorialproblem_sequential.hh + ../../tutorial/tutorial_sequential/tutorialspatialparams_sequential.hh) set(TEX_IMAGES PNG/box_disc.png diff --git a/dumux/common/intersectionmapper.hh b/dumux/common/intersectionmapper.hh index dcee8ab505f7253f58a363cab865c383a6ca760f..fadbd0422fd8400c7cf98665c2566a09b55777d7 100644 --- a/dumux/common/intersectionmapper.hh +++ b/dumux/common/intersectionmapper.hh @@ -24,6 +24,8 @@ #include #include +#include + /*! * \file * \brief defines an intersection mapper for mapping of global DOFs assigned @@ -44,7 +46,11 @@ class IntersectionMapper enum {dim=Grid::dimension}; typedef typename Grid::template Codim<0>::Entity Element; typedef typename GridView::Intersection Intersection; - typedef Dune::MultipleCodimMultipleGeomTypeMapper ElementMapper; +#if DUNE_VERSION_NEWER(DUNE_COMMON,2,6) + using ElementMapper = Dune::MultipleCodimMultipleGeomTypeMapper; +#else + using ElementMapper = Dune::MultipleCodimMultipleGeomTypeMapper; +#endif public: IntersectionMapper(const GridView& gridview) diff --git a/dumux/implicit/model.hh b/dumux/implicit/model.hh index 8edcd96150a1bb90522488e574c479ff46116af2..c8c531c91f20406c95bc899e2acebd57f65a47ab 100644 --- a/dumux/implicit/model.hh +++ b/dumux/implicit/model.hh @@ -31,6 +31,8 @@ #include #include +#include + namespace Dumux { @@ -69,7 +71,6 @@ class ImplicitModel typedef typename GridView::template Codim<0>::Entity Element; typedef typename Dune::ReferenceElements ReferenceElements; - typedef typename Dune::ReferenceElement ReferenceElement; enum { isBox = GET_PROP_VALUE(TypeTag, ImplicitIsBox) }; enum { dofCodim = isBox ? dim : 0 }; @@ -959,8 +960,12 @@ protected: for (const auto& element : elements(gridView_())) { Dune::GeometryType geomType = element.geometry().type(); - const ReferenceElement &refElement = ReferenceElements::general(geomType); +#if DUNE_VERSION_NEWER(DUNE_COMMON,2,6) + const auto refElement = ReferenceElements::general(geomType); +#else + const auto &refElement = ReferenceElements::general(geomType); +#endif for (const auto& intersection : intersections(gridView_(), element)) { if (intersection.boundary()) { if (isBox) diff --git a/dumux/implicit/problem.hh b/dumux/implicit/problem.hh index 81e9c326597abafd8cabf469588c1d76d38718fb..d940ba1e144178feaa4a36aa95e21cc45f752361 100644 --- a/dumux/implicit/problem.hh +++ b/dumux/implicit/problem.hh @@ -30,6 +30,8 @@ #include #include +#include + namespace Dumux { /*! @@ -103,11 +105,17 @@ public: : gridView_(gridView) , bBoxMin_(std::numeric_limits::max()) , bBoxMax_(-std::numeric_limits::max()) +#if DUNE_VERSION_NEWER(DUNE_COMMON,2,6) + , elementMapper_(gridView, Dune::mcmgElementLayout()) + , vertexMapper_(gridView, Dune::mcmgVertexLayout()) +#else , elementMapper_(gridView) , vertexMapper_(gridView) +#endif , timeManager_(&timeManager) , newtonMethod_(asImp_()) , newtonCtl_(asImp_()) + { // calculate the bounding box of the local partition of the grid view for (const auto& vertex : vertices(gridView)) { diff --git a/dumux/implicit/propertydefaults.hh b/dumux/implicit/propertydefaults.hh index 8cc5e846ac9b8deaaaed63d62d56cb42a41489f9..b83e553f49fdc37db8df41e111aa72bc88382989 100644 --- a/dumux/implicit/propertydefaults.hh +++ b/dumux/implicit/propertydefaults.hh @@ -41,6 +41,8 @@ #include "localjacobian.hh" #include "volumevariables.hh" +#include + namespace Dumux { // forward declarations @@ -71,16 +73,29 @@ SET_TYPE_PROP(ImplicitBase, NewtonMethod, NewtonMethod); SET_TYPE_PROP(ImplicitBase, NewtonController, NewtonController); //! Mapper for the grid view's vertices. +#if DUNE_VERSION_NEWER(DUNE_COMMON,2,6) +SET_TYPE_PROP(ImplicitBase, + VertexMapper, + Dune::MultipleCodimMultipleGeomTypeMapper); +#else SET_TYPE_PROP(ImplicitBase, VertexMapper, Dune::MultipleCodimMultipleGeomTypeMapper); +#endif + //! Mapper for the grid view's elements. +#if DUNE_VERSION_NEWER(DUNE_COMMON,2,6) +SET_TYPE_PROP(ImplicitBase, + ElementMapper, + Dune::MultipleCodimMultipleGeomTypeMapper); +#else SET_TYPE_PROP(ImplicitBase, ElementMapper, Dune::MultipleCodimMultipleGeomTypeMapper); +#endif //! Set the BaseModel to ImplicitModel SET_TYPE_PROP(ImplicitBase, BaseModel, ImplicitModel); diff --git a/dumux/io/gridcreator.hh b/dumux/io/gridcreator.hh index 94f18994a79f1454d027447f1a5ba25706595449..c681a5d29b55b7619889bf68a894577aaa1d5854 100644 --- a/dumux/io/gridcreator.hh +++ b/dumux/io/gridcreator.hh @@ -35,6 +35,7 @@ #include #include #include +#include #include #include #include @@ -63,8 +64,12 @@ // FoamGrid specific includes #if HAVE_DUNE_FOAMGRID #include +#if DUNE_VERSION_NEWER(DUNE_COMMON,2,6) +#include +#else #include #endif +#endif #include #include diff --git a/dumux/io/vtkmultiwriter.hh b/dumux/io/vtkmultiwriter.hh index cb7597643f129578d7e8fa2c0e49250915f45203..c7e67441953eca54848f0ada55ccd8fa869c5060 100644 --- a/dumux/io/vtkmultiwriter.hh +++ b/dumux/io/vtkmultiwriter.hh @@ -38,6 +38,8 @@ #include +#include + #if HAVE_MPI #include #endif @@ -56,8 +58,13 @@ class VtkMultiWriter { enum { dim = GridView::dimension }; +#if DUNE_VERSION_NEWER(DUNE_COMMON,2,6) + typedef Dune::MultipleCodimMultipleGeomTypeMapper VertexMapper; + typedef Dune::MultipleCodimMultipleGeomTypeMapper ElementMapper; +#else typedef Dune::MultipleCodimMultipleGeomTypeMapper VertexMapper; typedef Dune::MultipleCodimMultipleGeomTypeMapper ElementMapper; +#endif public: typedef Dune::VTKWriter VtkWriter; @@ -65,8 +72,13 @@ public: const std::string &simName = "", std::string multiFileName = "") : gridView_(gridView) +#if DUNE_VERSION_NEWER(DUNE_COMMON,2,6) + , elementMapper_(gridView, Dune::mcmgElementLayout()) + , vertexMapper_(gridView, Dune::mcmgVertexLayout()) +#else , elementMapper_(gridView) , vertexMapper_(gridView) +#endif { simName_ = (simName.empty())?"sim":simName; multiFileName_ = multiFileName; diff --git a/dumux/material/spatialparams/gstatrandomfield.hh b/dumux/material/spatialparams/gstatrandomfield.hh index c7dbe5e552a26776e65d54f4752e5235b58c2f7b..fbe0924f757563a7d072c36504f13d395019d3eb 100644 --- a/dumux/material/spatialparams/gstatrandomfield.hh +++ b/dumux/material/spatialparams/gstatrandomfield.hh @@ -30,6 +30,8 @@ #include #include +#include + namespace Dumux { @@ -54,7 +56,11 @@ class GstatRandomField using DataVector = std::vector; using Element = typename GridView::Traits::template Codim<0>::Entity; +#if DUNE_VERSION_NEWER(DUNE_COMMON,2,6) + using ElementMapper = Dune::MultipleCodimMultipleGeomTypeMapper; +#else using ElementMapper = Dune::MultipleCodimMultipleGeomTypeMapper; +#endif public: // Add field types if you want to implement e.g. tensor permeabilities. @@ -66,7 +72,12 @@ public: * \param gridView the used gridView */ GstatRandomField(const GridView& gridView) - : gridView_(gridView), elementMapper_(gridView), + : gridView_(gridView), +#if DUNE_VERSION_NEWER(DUNE_COMMON,2,6) + elementMapper_(gridView, Dune::mcmgElementLayout()), +#else + elementMapper_(gridView), +#endif data_(gridView.size(0)) {} /*! diff --git a/dumux/porousmediumflow/implicit/velocityoutput.hh b/dumux/porousmediumflow/implicit/velocityoutput.hh index 8f386f10425bf0e767c96c664a50cd7360ead310..91c7226105932640fffc43d5b36c365daea31366 100644 --- a/dumux/porousmediumflow/implicit/velocityoutput.hh +++ b/dumux/porousmediumflow/implicit/velocityoutput.hh @@ -31,6 +31,8 @@ #include #include +#include + namespace Dumux { @@ -66,7 +68,6 @@ class ImplicitVelocityOutput typedef Dune::FieldVector GlobalPosition; typedef typename Dune::ReferenceElements ReferenceElements; - typedef typename Dune::ReferenceElement ReferenceElement; enum { isBox = GET_PROP_VALUE(TypeTag, ImplicitIsBox) }; enum { dofCodim = isBox ? dim : 0 }; @@ -145,11 +146,14 @@ public: const auto geometry = element.geometry(); Dune::GeometryType geomType = geometry.type(); - const ReferenceElement &referenceElement - = ReferenceElements::general(geomType); +#if DUNE_VERSION_NEWER(DUNE_COMMON,2,6) + const auto refElement = ReferenceElements::general(geomType); +#else + const auto &refElement = ReferenceElements::general(geomType); +#endif const Dune::FieldVector& localPos - = referenceElement.position(0, 0); + = refElement.position(0, 0); // get the transposed Jacobian of the element mapping const typename Element::Geometry::JacobianTransposed jacobianT2 = diff --git a/dumux/porousmediumflow/sequential/properties.hh b/dumux/porousmediumflow/sequential/properties.hh index 09cfaa0a7abf0528e0a8df2f8a6a2a658a424a59..205cdaf30ca53f7efdb8d417a84743490939ce7c 100644 --- a/dumux/porousmediumflow/sequential/properties.hh +++ b/dumux/porousmediumflow/sequential/properties.hh @@ -24,6 +24,8 @@ #include #include +#include + /*! * \ingroup Sequential * \ingroup IMPETProperties @@ -136,12 +138,20 @@ public: /*! * \brief Mapper for the grid view's vertices. */ +#if DUNE_VERSION_NEWER(DUNE_COMMON,2,6) + using VertexMapper = Dune::MultipleCodimMultipleGeomTypeMapper; +#else typedef Dune::MultipleCodimMultipleGeomTypeMapper VertexMapper; +#endif /*! * \brief Mapper for the grid view's elements. */ +#if DUNE_VERSION_NEWER(DUNE_COMMON,2,6) + using ElementMapper = Dune::MultipleCodimMultipleGeomTypeMapper; +#else typedef Dune::MultipleCodimMultipleGeomTypeMapper ElementMapper; +#endif /*! * \brief The type of a solution at a fixed time. diff --git a/dumux/porousmediumflow/sequential/variableclass.hh b/dumux/porousmediumflow/sequential/variableclass.hh index 25a1e17554f3a892a13ccec99bf92a3d6b2cae59..289b1659fd1acb7b6d0d40c80747cdcbb1d9828b 100644 --- a/dumux/porousmediumflow/sequential/variableclass.hh +++ b/dumux/porousmediumflow/sequential/variableclass.hh @@ -21,6 +21,8 @@ #include "properties.hh" +#include + // for parallelization //#include @@ -77,10 +79,16 @@ public: * @param gridView a DUNE gridview object corresponding to diffusion and transport equation */ VariableClass(const GridView& gridView) : - gridView_(gridView), elementMapper_(gridView), vertexMapper_(gridView) + gridView_(gridView), +#if DUNE_VERSION_NEWER(DUNE_COMMON,2,6) + elementMapper_(gridView, Dune::mcmgElementLayout()), + vertexMapper_(gridView, Dune::mcmgVertexLayout()) +#else + elementMapper_(gridView), + vertexMapper_(gridView) +#endif {} - //! Initializes the variable class /*! Method initializes the cellData vector. * Should be called from problem init() diff --git a/test/geomechanics/el2p/el2pproblem.hh b/test/geomechanics/el2p/el2pproblem.hh index ae4f7fe3feca429a5fa025468a5be3215b39807a..f87399f2e411b405035e2678850471f5f6ffaa60 100644 --- a/test/geomechanics/el2p/el2pproblem.hh +++ b/test/geomechanics/el2p/el2pproblem.hh @@ -34,6 +34,8 @@ #include "el2pco2tables.hh" #include "el2pspatialparams.hh" +#include + namespace Dumux { template @@ -764,7 +766,14 @@ public: * * \param gridView The grid view */ - InitialPressSat(const GridView & gridView) : BaseT(gridView) , gridView_(gridView), vertexMapper_(gridView) + InitialPressSat(const GridView & gridView) + : BaseT(gridView) + , gridView_(gridView) +#if DUNE_VERSION_NEWER(DUNE_COMMON,2,6) + , vertexMapper_(gridView, Dune::mcmgVertexLayout()) +#else + , vertexMapper_(gridView) +#endif { // resize the pressure field vector with the number of vertices pInit_.resize(gridView.size(GridView::dimension)); diff --git a/tutorial/CMakeLists.txt b/tutorial/CMakeLists.txt index adc6d64473192a66d91670642a857fc5d54d8a88..72c0c32f694edd7ff50e4215e7f6f308fa696503 100644 --- a/tutorial/CMakeLists.txt +++ b/tutorial/CMakeLists.txt @@ -1,17 +1,5 @@ -add_input_file_links() - -add_dumux_test(tutorial_sequential tutorial_sequential tutorial_sequential.cc - ${CMAKE_CURRENT_BINARY_DIR}/tutorial_sequential) - -add_dumux_test(tutorial_implicit tutorial_implicit tutorial_implicit.cc - ${CMAKE_CURRENT_BINARY_DIR}/tutorial_implicit) - -#install sources -install(FILES -tutorial_implicit.cc -tutorial_sequential.cc -tutorialproblem_implicit.hh -tutorialproblem_sequential.hh -tutorialspatialparams_implicit.hh -tutorialspatialparams_sequential.hh -DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dumux/tutorial) +add_subdirectory(ex1) +add_subdirectory(ex2) +add_subdirectory(ex3) +add_subdirectory(tutorial_implicit) +add_subdirectory(tutorial_sequential) diff --git a/tutorial/ex1/CMakeLists.txt b/tutorial/ex1/CMakeLists.txt new file mode 100644 index 0000000000000000000000000000000000000000..8f7d6e6d612bb4367ac59616ff6b528303ed2a8b --- /dev/null +++ b/tutorial/ex1/CMakeLists.txt @@ -0,0 +1,10 @@ +# the immiscible two-phase simulation program +dune_add_test(NAME exercise1_2p + SOURCES exercise1_2p.cc) + +# the compositional two-phase simulation program +dune_add_test(NAME exercise1_2p2c + SOURCES exercise1_2p2c.cc) + +# add a symlink for the input file +dune_symlink_to_source_files(FILES "exercise1.input") diff --git a/tutorial/ex1/README.md b/tutorial/ex1/README.md new file mode 100644 index 0000000000000000000000000000000000000000..ca823871d5174ac78a3a5332cb3d1d1c612d0dfb --- /dev/null +++ b/tutorial/ex1/README.md @@ -0,0 +1,152 @@ +# Exercise #1 (DuMuX course) +
+## Problem set-up + +N2 is injected in an aquifer previously saturated with water with an injection rate of 0.001 kg/(s*m$`^2`$). +The aquifer is situated 2700 m below see level and the domain size is 60 m x 40 m. It consists of two layers, a moderately permeable one ($`\Omega_1`$) and a lower permeable one ($`\Omega_2`$). + + + +## Preparing the exercise + +* Navigate to the directory `dumux/tutorial/ex1` + +_Exercise 1_ deals with two problems: a two-phase immiscible problem (__2p__) and a two-phase compositional problem (__2p2c__). They both set up the same scenario with the difference that the 2p2c assumes a miscible fluid state for the two fluids (water and gaseous N2) and the 2p model assumes an immiscible fluid state. + +

+### Task 1: Getting familiar with the code +
+ +Locate all the files you will need for this exercise +* The __main file__ for the __2p__ problem: `exercise1_2p.cc` +* The __problem file__ for the __2p__ problem: `injection2pproblem.hh` +* The __main file__ for the __2p2c__ problem: `exercise1_2p2c.cc` +* The problem file for the __2p2c__ problem: `injection2p2cproblem.hh` +* The shared __spatial parameters file__: `injection2pspatialparams.hh` +* The shared __input file__: `exercise1.input` + +


+### Task 2: Compiling and running an executable +
+ +* Change to the build-directory + +```bash +cd ../../build-cmake/tutorial/exercise1 +``` + +* Compile both executables `exercise1_2p` and `exercise1_2p2c` + +```bash +make exercise1_2p exercise1_2p2c +``` + +* Execute the two problems and inspect the result + +```bash +./exercise1_2p exercise1.input +./exercise1_2p2c exercise1.input +``` + +* you can look at the results with paraview + +```bash +paraview injection-2p2c.pvd +``` + +


+### Task 3: Changing input parameters +
+ +In the input file `exercise1.input` you can find the following section + +```ini +[SpatialParams] +EntryPressureFine = 4.5e4 +EntryPressureCoarse = 1e4 +``` + +* Change the values for the fine entry pressure in the input file to a lower value and compare the results with the previous solution. You don't need to recompile the executable. + +


+### Task 4: Runtime parameters +
+ +The injection rate is currently hard-coded in `injection2p2cproblem.hh` to $`1e-4 kg/(s m^2)`$. + +```c++ + // set the Neumann values for the Nitrogen component balance + // convert from units kg/(s*m^2) to mole/(s*m^2) +values[Indices::contiNEqIdx] = -1e-4/FluidSystem::molarMass(FluidSystem::nCompIdx); +values[Indices::contiWEqIdx] = 0.0; +``` + +We want to be able to set it at runtime. To this end, +* use the following DuMuX macro to read a runtime parameter from the input file + +```c++ +// read the injection rate from the input file at run time +const auto totalAreaSpecificInflow_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, , , ); +``` + +* Replace +``,``,`` by what is appropriate for your case: + * `` is the type of the parameter to read + * `` is the group in the input file + * `` is the name of the parameter in the input file + +Note that due to the way the macro works, the names are specified as plain text without "quotation marks". Follow the instructions given as a + +```c++ +// TODO: dumux-course-task +``` +in the `injection2p2cproblem.hh` file and also remember to also set the parameter totalAreaSpecificInflow in the input file. + +* Check the influence of that parameter on the simulation result by rerunning the simulation with different injection rates. Remember to also set the parameter totalAreaSpecificInflow in the input file. + +Again, you don't need to recompile the program. + +


+### 5. Setting up a new executable (for a non-isothermal simulation) +
+ +* Set up a new cc file called `exercise1_2pni.cc` by copying and renaming `exercise1_2p.cc` + +```bash +cp exercise1_2p.cc exercise1_2pni.cc +``` + +* In `exercise1_2pni.cc`, include the header `injection2pniproblem.hh` instead of the isothermal problem file `injection2pproblem.hh`. +* Add a new executable in `CMakeLists.txt` by adding the lines + +```cmake +dune_add_test(NAME injection2pniproblem + SOURCES injection2pniproblem.cc) +``` + +* Test that everything compiles without error + +```bash +make # should rerun cmake +make injection2pniproblem # builds new executable +``` + +


+### 6. Setting up a non-isothermal __2pni__ test problem +
+ +* Open the file `injection2pniproblem.hh`. It is a copy of the `injection2pproblem.hh` with some useful comments on how to implement a non-isothermal model. Look for comments containing + +```c++ +// TODO: dumux-course-task +``` + +* The following set-up should be realized: + + __Boundary conditions:__ Dirichlet conditions for the temperature with a temperature gradient of 0.03 K/m and a starting temperature of 283 K. + + __Initial conditions:__ The same temperature gradient as in the boundary conditions with an additional lens (with position: 20 < x < 30, 5 < y < 35), which has an initial temperature of 380 K. + + + +The non-isothermal model requires additional spatial parameters like the thermal conductivity. They are already implemented in `injection2pspatialparams.hh`, you just need to _uncomment_ them. diff --git a/tutorial/ex1/exercise1.input b/tutorial/ex1/exercise1.input new file mode 100644 index 0000000000000000000000000000000000000000..4f9a57e69a8ccd12d4831c321115806d7c7a92f5 --- /dev/null +++ b/tutorial/ex1/exercise1.input @@ -0,0 +1,24 @@ +[TimeManager] +DtInitial = 3600 # in seconds +TEnd = 3.154e9 # in seconds, i.e ten years + +[Grid] +LowerLeft = 0 0 +UpperRight = 60 40 +Cells = 24 16 + +[Problem] +Name = injection +OnlyPlotMaterialLaws = true +AquiferDepth = 2700.0 # m +InjectionDuration = 2.628e6 # in seconds, i.e. one month + +#TODO: dumux-course-task: +#set totalAreaSpecificInflow + +[SpatialParams] +PermeabilityAquitard = 1e-15 # m^2 +EntryPressureAquitard = 4.5e4 # Pa +PermeabilityAquifer = 1e-12 # m^2 +EntryPressureAquifer = 1e4 # Pa + diff --git a/tutorial/ex1/exercise1_2p.cc b/tutorial/ex1/exercise1_2p.cc new file mode 100644 index 0000000000000000000000000000000000000000..6152aa5d204c851134b078f0d4e46fc030353d50 --- /dev/null +++ b/tutorial/ex1/exercise1_2p.cc @@ -0,0 +1,64 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/***************************************************************************** + * See the file COPYING for full copying permissions. * + * * + * This program is free software: you can redistribute it and/or modify * + * it under the terms of the GNU General Public License as published by * + * the Free Software Foundation, either version 2 of the License, or * + * (at your option) any later version. * + * * + * This program is distributed in the hope that it will be useful, * + * but WITHOUT ANY WARRANTY; without even the implied warranty of * + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * + * GNU General Public License for more details. * + * * + * You should have received a copy of the GNU General Public License * + * along with this program. If not, see . * + *****************************************************************************/ +/*! + * \file + * + * \brief DOC ME! + */ +#include "config.h" +#include "injection2pproblem.hh" +#include + +/*! + * \brief Provides an interface for customizing error messages associated with + * reading in parameters. + * + * \param progName The name of the program, that was tried to be started. + * \param errorMsg The error message that was issued by the start function. + * Comprises the thing that went wrong and a general help message. + */ +void usage(const char *progName, const std::string &errorMsg) +{ + if (errorMsg.size() > 0) { + std::string errorMessageOut = "\nUsage: "; + errorMessageOut += progName; + errorMessageOut += " [options]\n"; + errorMessageOut += errorMsg; + errorMessageOut += "\n\nThe list of mandatory options for this program is:\n" + "\t-TimeManager.TEnd End of the simulation [s] \n" + "\t-TimeManager.DtInitial Initial timestep size [s] \n" + "\t-Grid.File Name of the file containing the grid \n" + "\t definition in DGF format\n" + "\t-SpatialParams.LensLowerLeft coordinates of the lower left corner of the lens [m] \n" + "\t-SpatialParams.LensUpperRight coordinates of the upper right corner of the lens [m] \n" + "\t-Problem.Name String for naming of the output files \n" + "\n"; + std::cout << errorMessageOut + << "\n"; + } +} + +//////////////////////// +// the main function +//////////////////////// +int main(int argc, char** argv) +{ + typedef TTAG(InjectionCCProblem2P) ProblemTypeTag; + return Dumux::start(argc, argv, usage); +} diff --git a/tutorial/ex1/exercise1_2p2c.cc b/tutorial/ex1/exercise1_2p2c.cc new file mode 100644 index 0000000000000000000000000000000000000000..aba0b0e8c76986dd481d8ad9978ddb7874cb01ae --- /dev/null +++ b/tutorial/ex1/exercise1_2p2c.cc @@ -0,0 +1,67 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/***************************************************************************** + * See the file COPYING for full copying permissions. * + * * + * This program is free software: you can redistribute it and/or modify * + * it under the terms of the GNU General Public License as published by * + * the Free Software Foundation, either version 2 of the License, or * + * (at your option) any later version. * + * * + * This program is distributed in the hope that it will be useful, * + * but WITHOUT ANY WARRANTY; without even the implied warranty of * + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * + * GNU General Public License for more details. * + * * + * You should have received a copy of the GNU General Public License * + * along with this program. If not, see . * + *****************************************************************************/ +/*! + * \file + * + * \brief Test for the two-phase two-component CC model. + */ +#include +#include "injection2p2cproblem.hh" +#include + +/*! + * \brief Provides an interface for customizing error messages associated with + * reading in parameters. + * + * \param progName The name of the program, that was tried to be started. + * \param errorMsg The error message that was issued by the start function. + * Comprises the thing that went wrong and a general help message. + */ +void usage(const char *progName, const std::string &errorMsg) +{ + if (errorMsg.size() > 0) { + std::string errorMessageOut = "\nUsage: "; + errorMessageOut += progName; + errorMessageOut += " [options]\n"; + errorMessageOut += errorMsg; + errorMessageOut += "\n\nThe list of mandatory options for this program is:\n" + "\t-TimeManager.TEnd End of the simulation [s] \n" + "\t-TimeManager.DtInitial Initial timestep size [s] \n" + "\t-Grid.File Name of the file containing the grid \n" + "\t definition in DGF format\n" + "\t-FluidSystem.NTemperature Number of tabularization entries [-] \n" + "\t-FluidSystem.NPressure Number of tabularization entries [-] \n" + "\t-FluidSystem.PressureLow Low end for tabularization of fluid properties [Pa] \n" + "\t-FluidSystem.PressureHigh High end for tabularization of fluid properties [Pa] \n" + "\t-FluidSystem.TemperatureLow Low end for tabularization of fluid properties [Pa] \n" + "\t-FluidSystem.TemperatureHigh High end for tabularization of fluid properties [Pa] \n" + "\t-SimulationControl.Name The name of the output files [-] \n" + "\t-InitialConditions.Temperature Initial temperature in the reservoir [K] \n" + "\t-InitialConditions.DepthBOR Depth below ground surface [m] \n"; + + std::cout << errorMessageOut + << "\n"; + } +} + +int main(int argc, char** argv) +{ + typedef TTAG(Injection2p2pcCCProblem) ProblemTypeTag; + return Dumux::start(argc, argv, usage); +} diff --git a/tutorial/ex1/injection2p2cproblem.hh b/tutorial/ex1/injection2p2cproblem.hh new file mode 100644 index 0000000000000000000000000000000000000000..aada803f77f5bcfbcb06a4e61122de6a7a2c560a --- /dev/null +++ b/tutorial/ex1/injection2p2cproblem.hh @@ -0,0 +1,307 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/***************************************************************************** + * See the file COPYING for full copying permissions. * + * * + * This program is free software: you can redistribute it and/or modify * + * it under the terms of the GNU General Public License as published by * + * the Free Software Foundation, either version 2 of the License, or * + * (at your option) any later version. * + * * + * This program is distributed in the hope that it will be useful, * + * but WITHOUT ANY WARRANTY; without even the implied warranty of * + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * + * GNU General Public License for more details. * + * * + * You should have received a copy of the GNU General Public License * + * along with this program. If not, see . * + *****************************************************************************/ +/*! + * \file + * + * \brief Problem where air is injected under a low permeable layer in a depth of 2700m. + */ +#ifndef DUMUX_INJECTION_2P2C_PROBLEM_HH +#define DUMUX_INJECTION_2P2C_PROBLEM_HH + +#include +#include +#include + +#include "injection2pspatialparams.hh" + +namespace Dumux +{ + +template +class Injection2p2cProblem; + +namespace Properties +{ +NEW_TYPE_TAG(Injection2p2cProblem, INHERITS_FROM(TwoPTwoC, InjectionSpatialParams)); +NEW_TYPE_TAG(Injection2p2cBoxProblem, INHERITS_FROM(BoxModel, Injection2p2cProblem)); +NEW_TYPE_TAG(Injection2p2pcCCProblem, INHERITS_FROM(CCModel, Injection2p2cProblem)); + +// Set the grid type +SET_TYPE_PROP(Injection2p2cProblem, Grid, Dune::YaspGrid<2>); + +// Set the problem property +SET_TYPE_PROP(Injection2p2cProblem, Problem, Injection2p2cProblem); + +// Set fluid configuration +SET_TYPE_PROP(Injection2p2cProblem, + FluidSystem, + FluidSystems::H2ON2); + +// Define whether mole(true) or mass (false) fractions are used +SET_BOOL_PROP(Injection2p2cProblem, UseMoles, true); +} + + +/*! + * \ingroup TwoPTwoCModel + * \ingroup ImplicitTestProblems + * \brief Problem where air is injected under a low permeable layer in a depth of 2700m. + * + * The domain is sized 60m times 40m and consists of two layers, a moderately + * permeable one (\f$ K=10e-12\f$) and one with a lower permeablility (\f$ K=10e-13\f$) + * in the rest of the domain. + * + * Nitrogen is injected into a water-filled aquifer through a well. First, we inject for one month. + * Then, we continue simulating the development of the nitrogen plume for 10 years.This is realized with a solution-dependent Neumann boundary condition at the right boundary + * (\f$ 7m./exercise1_2p2c -ParameterFile exercise1.input> + */ +template +class Injection2p2cProblem : public ImplicitPorousMediaProblem +{ + typedef ImplicitPorousMediaProblem ParentType; + + typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView; + typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar; + typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem; + + // grid world dimension + static constexpr auto dimWorld = GridView::dimensionworld; + + typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices; + typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables; + typedef typename GET_PROP_TYPE(TypeTag, BoundaryTypes) BoundaryTypes; + typedef typename GET_PROP_TYPE(TypeTag, TimeManager) TimeManager; + typedef typename GET_PROP_TYPE(TypeTag, MaterialLaw) MaterialLaw; + typedef Dune::FieldVector GlobalPosition; + +public: + /*! + * \brief The constructor + * + * \param timeManager The time manager + * \param gridView The grid view + */ + Injection2p2cProblem(TimeManager &timeManager, + const GridView &gridView) + : ParentType(timeManager, gridView) + { + // initialize the tables of the fluid system + FluidSystem::init(/*tempMin=*/273.15, + /*tempMax=*/423.15, + /*numTemp=*/50, + /*pMin=*/0.0, + /*pMax=*/30e6, + /*numP=*/300); + + // name of the problem and output file + name_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, std::string, Problem, Name); + // depth of the aquifer, units: m + aquiferDepth_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Problem, AquiferDepth); + // the duration of the injection, units: second + injectionDuration_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Problem, InjectionDuration); + + // TODO: dumux-course-task + // Get the specific inflow of 1e-4 kg/(s m^2) from the input file (totalAreaSpecificInflow_) here as it is done for the injectionDuration_. + + } + + /*! + * \brief User defined output after the time integration + * + * Will be called diretly after the time integration. + */ + void postTimeStep() + { + // Calculate storage terms + PrimaryVariables storageW, storageN; + this->model().globalPhaseStorage(storageW, Indices::wPhaseIdx); + this->model().globalPhaseStorage(storageN, Indices::nPhaseIdx); + + // Write mass balance information for rank 0 + if (this->gridView().comm().rank() == 0) { + std::cout<<"Storage: wetting=[" << storageW << "]" + << " nonwetting=[" << storageN << "]\n"; + } + } + + + /*! + * \name Problem parameters + */ + // \{ + + /*! + * \brief Returns the problem name + * + * This is used as a prefix for files generated by the simulation. + */ + std::string name() const + { return name_+"-2p2c"; } + + /*! + * \brief Returns the temperature \f$ K \f$ + */ + Scalar temperature() const + { + return 273.15 + 30; // [K] + } + + /*! + * \brief Returns the source term + * + * \param values Stores the source values for the conservation equations in + * \f$ [ \textnormal{unit of primary variable} / (m^\textrm{dim} \cdot s )] \f$ + * \param globalPos The global position + */ + void sourceAtPos(PrimaryVariables &values, + const GlobalPosition &globalPos) const + { + values = 0; + } + + // \} + + /*! + * \name Boundary conditions + */ + // \{ + + /*! + * \brief Specifies which kind of boundary condition should be + * used for which equation on a given boundary segment + * + * \param values Stores the value of the boundary type + * \param globalPos The global position + */ + void boundaryTypesAtPos(BoundaryTypes &values, + const GlobalPosition &globalPos) const + { + if (globalPos[dimWorld-1] < eps_) + values.setAllDirichlet(); + + // and Neuman boundary conditions everywhere else + // note that we don't differentiate between Neumann and Robin boundary types + else + values.setAllNeumann(); + } + + /*! + * \brief Evaluates the boundary conditions for a Dirichlet + * boundary segment + * + * \param values Stores the Dirichlet values for the conservation equations in + * \f$ [ \textnormal{unit of primary variable} ] \f$ + * \param globalPos The global position + */ + void dirichletAtPos(PrimaryVariables &values, const GlobalPosition &globalPos) const + { + initialAtPos(values, globalPos); + } + + /*! + * \brief Evaluates the boundary conditions for a Neumann boundary segment. + * + * \param values Stores the Neumann values for the conservation equations in + * \f$ [ \textnormal{unit of conserved quantity} / (m^(dim-1) \cdot s )] \f$ + * \param globalPos The globalPosition of the boundary interface + */ + void neumannAtPos(PrimaryVariables &values, + const GlobalPosition &globalPos) const + { + // initialize values to zero, i.e. no-flow Neumann boundary conditions + values = 0; + + //if we are inside the injection zone set inflow Neumann boundary conditions + if (this->timeManager().time() + this->timeManager().timeStepSize() < injectionDuration_ + && globalPos[1] < 15 + eps_ && globalPos[1] > 7 - eps_ && globalPos[0] > 0.9*this->bBoxMax()[0]) + { + // TODO: dumux-course-task + //instead of setting -1e-4 here directly use totalAreaSpecificInflow_ in the computation + + // set the Neumann values for the Nitrogen component balance + // convert from units kg/(s*m^2) to mole/(s*m^2) + values[Indices::contiNEqIdx] = -1e-4/FluidSystem::molarMass(FluidSystem::nCompIdx); + values[Indices::contiWEqIdx] = 0.0; + } + } + + // \} + + /*! + * \name Volume terms + */ + // \{ + + /*! + * \brief Evaluates the initial values for a control volume + * + * \param values Stores the initial values for the conservation equations in + * \f$ [ \textnormal{unit of primary variables} ] \f$ + * \param globalPos The global position + */ + void initialAtPos(PrimaryVariables &values, const GlobalPosition &globalPos) const + { + // get the water density at atmospheric conditions + const Scalar densityW = FluidSystem::H2O::liquidDensity(temperature(), 1.0e5); + + // assume an intially hydrostatic liquid pressure profile + // note: we subtract rho_w*g*h because g is defined negative + const Scalar pw = 1.0e5 - densityW*this->gravity()[dimWorld-1]*(aquiferDepth_ - globalPos[dimWorld-1]); + + // initially we have some nitrogen dissolved + // saturation mole fraction would be + // moleFracLiquidN2 = (pw + pc + p_vap^sat)/henry; + const Scalar moleFracLiquidN2 = pw*0.95/BinaryCoeff::H2O_N2::henry(temperature()); + + // note that because we start with a single phase system the primary variables + // are pl and x^w_N2. This will switch as soon after we start injecting to a two + // phase system so the primary variables will be pl and Sn (non-wetting saturation). + values[Indices::switchIdx] = moleFracLiquidN2; + values[Indices::pressureIdx] = pw; + } + + /*! + * \brief Return the initial phase state inside a control volume. + * + * \param globalPos The global position + */ + int initialPhasePresenceAtPos(const GlobalPosition &globalPos) const + { return Indices::wPhaseOnly; } + + // \} + +private: + static constexpr Scalar eps_ = 1e-6; + std::string name_; //! Problem name + Scalar aquiferDepth_; //! Depth of the aquifer in m + Scalar injectionDuration_; //! Duration of the injection in seconds + //TODO: dumux-course-task + //define the Scalar totalAreaSpecificInflow_ here + +}; +} //end namespace + +#endif diff --git a/tutorial/ex1/injection2pniproblem.hh b/tutorial/ex1/injection2pniproblem.hh new file mode 100644 index 0000000000000000000000000000000000000000..f7453851848625ae1af19337cac0e930887efd4d --- /dev/null +++ b/tutorial/ex1/injection2pniproblem.hh @@ -0,0 +1,304 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/***************************************************************************** + * See the file COPYING for full copying permissions. * + * * + * This program is free software: you can redistribute it and/or modify * + * it under the terms of the GNU General Public License as published by * + * the Free Software Foundation, either version 2 of the License, or * + * (at your option) any later version. * + * * + * This program is distributed in the hope that it will be useful, * + * but WITHOUT ANY WARRANTY; without even the implied warranty of * + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * + * GNU General Public License for more details. * + * * + * You should have received a copy of the GNU General Public License * + * along with this program. If not, see . * + *****************************************************************************/ +/*! + * \file + * + * \brief Non-isothermal gas injection problem where a gas (e.g. air) is injected into a fully + * water saturated medium. During buoyancy driven upward migration the gas + * passes a high temperature area. + */ + +#ifndef DUMUX_INJECTION_PROBLEM_2PNI_HH +#define DUMUX_INJECTION_PROBLEM_2PNI_HH + +#include +#include + +#include + +#include "injection2pspatialparams.hh" + +namespace Dumux { + +template +class InjectionProblem2PNI; + +namespace Properties +{ + /*! +* TODO:dumux-course-task: +* inherit from the TwoPNI model instead of TwoP here +*/ +NEW_TYPE_TAG(InjectionProblem2PNI, INHERITS_FROM(TwoP, InjectionSpatialParams)); +NEW_TYPE_TAG(InjectionCCProblem2PNI, INHERITS_FROM(CCModel, InjectionProblem2PNI)); + +// Set the grid type +SET_TYPE_PROP(InjectionProblem2PNI, Grid, Dune::YaspGrid<2>); + +// Set the problem property +SET_TYPE_PROP(InjectionProblem2PNI, Problem, InjectionProblem2PNI); + +// Use the same fluid system as the 2p2c injection problem +SET_TYPE_PROP(InjectionProblem2PNI, FluidSystem, FluidSystems::H2ON2); + +} + +/*! + * \ingroup TwoPModel + * \ingroup ImplicitTestProblems + * \brief Non-isothermal gas injection problem where a gas (e.g. air) is injected into a fully + * water saturated medium. During buoyancy driven upward migration the gas + * passes a high temperature area. + * + * The domain is sized 60 m times 40 m. The rectangular area with the increased temperature (380 K) + * starts at (20 m, 5 m) and ends at (30 m, 35 m) + * + * For the mass conservation equation neumann boundary conditions are used on + * the top, on the bottom and on the right of the domain, while dirichlet conditions + * apply on the left boundary. + * For the energy conservation equation dirichlet boundary conditions are applied + * on all boundaries. + * + * Gas is injected at the right boundary from 7 m to 15 m at a rate of + * 0.001 kg/(s m), the remaining neumann boundaries are no-flow + * boundaries. + * + * At the dirichlet boundaries a hydrostatic pressure, a gas saturation of zero and + * a geothermal temperature gradient of 0.03 K/m are applied. + * + * This problem uses the \ref TwoPModel and \ref NIModel model. + * + * To run the simulation execute the following line in shell: + * ./exercise1_2pni -ParameterFile exercise1.input + */ +template +class InjectionProblem2PNI : public ImplicitPorousMediaProblem +{ + typedef ImplicitPorousMediaProblem ParentType; + typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView; + typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar; + + + typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices; + enum { + pressureIdx = Indices::pressureIdx, + saturationIdx = Indices::saturationIdx, + + contiNEqIdx = Indices::contiNEqIdx, + + /* + * TODO:dumux-course-task: + * get the temperatureIdx as the index for the primary variable temperature and the + * energyIdx as the index for the energy conservation equation for your convinience. + */ + + // world dimension + dimWorld = GridView::dimensionworld + }; + + + typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables; + typedef typename GET_PROP_TYPE(TypeTag, BoundaryTypes) BoundaryTypes; + typedef typename GET_PROP_TYPE(TypeTag, TimeManager) TimeManager; + + typedef typename GridView::template Codim<0>::Entity Element; + typedef typename GridView::Intersection Intersection; + + typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry; + typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem; + + typedef Dune::FieldVector GlobalPosition; + + enum { isBox = GET_PROP_VALUE(TypeTag, ImplicitIsBox) }; + +public: + /*! + * \brief The constructor + * + * \param timeManager The time manager + * \param gridView The grid view + */ + InjectionProblem2PNI(TimeManager &timeManager, const GridView &gridView) + : ParentType(timeManager, gridView) + { + // initialize the tables of the fluid system + FluidSystem::init(/*tempMin=*/273.15, + /*tempMax=*/423.15, + /*numTemp=*/50, + /*pMin=*/0.0, + /*pMax=*/30e6, + /*numP=*/300); + + // name of the problem and output file + name_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, std::string, Problem, Name); + // depth of the aquifer, units: m + aquiferDepth_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Problem, AquiferDepth); + // inflow rate of nitrogen water vapor mixture, units: kg/(s m^2) + // the duration of the injection, units: second + injectionDuration_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Problem, InjectionDuration); + } + + /*! + * \name Problem parameters + */ + // \{ + + /*! + * \brief Returns the problem name + * + * This is used as a prefix for files generated by the simulation. + */ + const std::string name() const + { return "injection-2pni";} + + /*! + * \brief Returns the source term + * + * \param values Stores the source values for the conservation equations in + * \f$ [ \textnormal{unit of primary variable} / (m^\textrm{dim} \cdot s )] \f$ + * \param globalPos The global position + */ + void sourceAtPos(PrimaryVariables &values, + const GlobalPosition &globalPos) const + { + values = 0; + } + + // \} + + /*! + * \name Boundary conditions + */ + // \{ + + /*! + * \brief Specifies which kind of boundary condition should be + * used for which equation on a given boundary segment + * + * \param values Stores the value of the boundary type + * \param globalPos The global position + */ + void boundaryTypesAtPos(BoundaryTypes &values, + const GlobalPosition &globalPos) const + { + if (globalPos[0] < eps_) + values.setAllDirichlet(); + else + values.setAllNeumann(); + + /*! + * dumux-course-task: + * set dirichlet conditions for the temperature index everywhere + */ + } + + /*! + * \brief Evaluates the boundary conditions for a Dirichlet + * boundary segment + * + * \param values Stores the Dirichlet values for the conservation equations in + * \f$ [ \textnormal{unit of primary variable} ] \f$ + * \param globalPos The global position + */ + void dirichletAtPos(PrimaryVariables &values, const GlobalPosition &globalPos) const + { + // assume a constant density + const Scalar densityW = 1000; + + // assume an intially hydrostatic liquid pressure profile + // note: we subtract rho_w*g*h because g is defined negative + const Scalar pw = 1.0e5 - densityW*this->gravity()[dimWorld-1]*(aquiferDepth_ - globalPos[dimWorld-1]); + + values[Indices::pressureIdx] = pw; + values[saturationIdx] = 0.0; + /*! + * TODO:dumux-course-task: + * set a temperature gradient of 0.03 K per m beginning at 283 K here. + * Hint: you can use maxDepth_ and the globalPos similar to the pressure gradient + */ + } + + /*! + * \brief Evaluates the boundary conditions for a Neumann boundary segment. + * + * \param values Stores the Neumann values for the conservation equations in + * \f$ [ \textnormal{unit of conserved quantity} / (m^(dim-1) \cdot s )] \f$ + * \param globalPos The globalPosition of the boundary interface + */ + void neumannAtPos(PrimaryVariables &values, + const GlobalPosition &globalPos) const + { + // initialize values to zero, i.e. no-flow Neumann boundary conditions + values = 0; + values = 0; + + if (this->timeManager().time() + this->timeManager().timeStepSize() < injectionDuration_ + && globalPos[1] < 15 + eps_ && globalPos[1] > 7 - eps_ && globalPos[0] > 0.9*this->bBoxMax()[0]) + { + // inject air. negative values mean injection + values[Indices::contiNEqIdx] = -1e-4; // kg/(s*m^2) + values[Indices::contiWEqIdx] = 0.0; + } + } + + // \} + + + /*! + * \name Volume terms + */ + // \{ + + /*! + * \brief Evaluates the initial values for a control volume + * + * \param values Stores the initial values for the conservation equations in + * \f$ [ \textnormal{unit of primary variables} ] \f$ + * \param globalPos The global position + */ + void initialAtPos(PrimaryVariables &values, const GlobalPosition &globalPos) const + { + // assume a constant density + const Scalar densityW = 1000; + + // assume an intially hydrostatic liquid pressure profile + // note: we subtract rho_w*g*h because g is defined negative + const Scalar pw = 1.0e5 - densityW*this->gravity()[dimWorld-1]*(aquiferDepth_ - globalPos[dimWorld-1]); + + values[Indices::pressureIdx] = pw; + values[saturationIdx] = 0.0; + + /*! + * TODO:dumux-course-task: + * set a temperature gradient of 0.03 K per m beginning at 283 K here. + * Hint: you can use maxDepth_ and the globalPos similar to the pressure gradient + * use globalPos[0] and globalpos[1] to implement the high temperature lens with 380 K + */ + } + + +private: + static constexpr Scalar eps_ = 1e-6; + std::string name_; //! Problem name + Scalar aquiferDepth_; //! Depth of the aquifer in m + Scalar injectionDuration_; //! Duration of the injection in seconds +}; +} //end namespace + +#endif diff --git a/tutorial/ex1/injection2pproblem.hh b/tutorial/ex1/injection2pproblem.hh new file mode 100644 index 0000000000000000000000000000000000000000..c990ffb09b4c2b41aa55834bda0c9a2c7180f0ef --- /dev/null +++ b/tutorial/ex1/injection2pproblem.hh @@ -0,0 +1,276 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/***************************************************************************** + * See the file COPYING for full copying permissions. * + * * + * This program is free software: you can redistribute it and/or modify * + * it under the terms of the GNU General Public License as published by * + * the Free Software Foundation, either version 2 of the License, or * + * (at your option) any later version. * + * * + * This program is distributed in the hope that it will be useful, * + * but WITHOUT ANY WARRANTY; without even the implied warranty of * + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * + * GNU General Public License for more details. * + * * + * You should have received a copy of the GNU General Public License * + * along with this program. If not, see . * + *****************************************************************************/ +/*! + * \file + * + * \brief Non-isothermal gas injection problem where a gas (e.g. air) is injected into a fully + * water saturated medium. During buoyancy driven upward migration the gas + * passes a high temperature area. + */ + +#ifndef DUMUX_INJECTION_PROBLEM_2P_HH +#define DUMUX_INJECTION_PROBLEM_2P_HH + +#include +#include + +#include + +#include "injection2pspatialparams.hh" + +namespace Dumux { + +template +class InjectionProblem2P; + +namespace Properties +{ +NEW_TYPE_TAG(InjectionProblem2P, INHERITS_FROM(TwoP, InjectionSpatialParams)); +NEW_TYPE_TAG(InjectionCCProblem2P, INHERITS_FROM(CCModel, InjectionProblem2P)); + +// Set the grid type +SET_TYPE_PROP(InjectionProblem2P, Grid, Dune::YaspGrid<2>); + +// Set the problem property +SET_TYPE_PROP(InjectionProblem2P, Problem, InjectionProblem2P); + +// Use the same fluid system as the 2p2c injection problem +SET_TYPE_PROP(InjectionProblem2P, FluidSystem, FluidSystems::H2ON2); + +} + +/*! + * \ingroup TwoPModel + * \ingroup ImplicitTestProblems + * \brief Gas injection problem where a gas (here nitrogen) is injected into a fully + * water saturated medium. During buoyancy driven upward migration the gas + * passes a high temperature area. + * + * The domain is sized 60 m times 40 m. + * + * For the mass conservation equation neumann boundary conditions are used on + * the top, on the bottom and on the right of the domain, while dirichlet conditions + * apply on the left boundary. + * + * Gas is injected at the right boundary from 7 m to 15 m at a rate of + * 0.001 kg/(s m), the remaining neumann boundaries are no-flow + * boundaries. + * + * At the dirichlet boundaries a hydrostatic pressure and a gas saturation of zero a + * + * This problem uses the \ref TwoPModel model. + * + * To run the simulation execute the following line in shell: + * ./exercise1_2p -ParameterFile exercise1.input + */ +template +class InjectionProblem2P : public ImplicitPorousMediaProblem +{ + typedef ImplicitPorousMediaProblem ParentType; + typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView; + typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar; + + + typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices; + enum { + pressureIdx = Indices::pressureIdx, + saturationIdx = Indices::saturationIdx, + + contiNEqIdx = Indices::contiNEqIdx, + + // world dimension + dimWorld = GridView::dimensionworld + }; + + + typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables; + typedef typename GET_PROP_TYPE(TypeTag, BoundaryTypes) BoundaryTypes; + typedef typename GET_PROP_TYPE(TypeTag, TimeManager) TimeManager; + + typedef typename GridView::template Codim<0>::Entity Element; + typedef typename GridView::Intersection Intersection; + + typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry; + typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem; + + typedef Dune::FieldVector GlobalPosition; + + enum { isBox = GET_PROP_VALUE(TypeTag, ImplicitIsBox) }; + +public: + /*! + * \brief The constructor + * + * \param timeManager The time manager + * \param gridView The grid view + */ + InjectionProblem2P(TimeManager &timeManager, const GridView &gridView) + : ParentType(timeManager, gridView) + { + // initialize the tables of the fluid system + FluidSystem::init(/*tempMin=*/273.15, + /*tempMax=*/423.15, + /*numTemp=*/50, + /*pMin=*/0.0, + /*pMax=*/30e6, + /*numP=*/300); + + // name of the problem and output file + name_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, std::string, Problem, Name); + // depth of the aquifer, units: m + aquiferDepth_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Problem, AquiferDepth); + // inflow rate of nitrogen water vapor mixture, units: kg/(s m^2) + // the duration of the injection, units: second + injectionDuration_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Problem, InjectionDuration); + } + + /*! + * \name Problem parameters + */ + // \{ + + /*! + * \brief Returns the problem name + * + * This is used as a prefix for files generated by the simulation. + */ + const std::string name() const + { return "injection-2p"; } + + /*! + * \brief Returns the temperature \f$ K \f$ + */ + Scalar temperature() const + { + return 273.15 + 30; // [K] + } + + + + /*! + * \brief Returns the source term + * + * \param values Stores the source values for the conservation equations in + * \f$ [ \textnormal{unit of primary variable} / (m^\textrm{dim} \cdot s )] \f$ + * \param globalPos The global position + */ + void sourceAtPos(PrimaryVariables &values, + const GlobalPosition &globalPos) const + { + values = 0; + } + + // \} + + /*! + * \name Boundary conditions + */ + // \{ + + /*! + * \brief Specifies which kind of boundary condition should be + * used for which equation on a given boundary segment + * + * \param values Stores the value of the boundary type + * \param globalPos The global position + */ + void boundaryTypesAtPos(BoundaryTypes &values, + const GlobalPosition &globalPos) const + { + if (globalPos[0] < eps_) + values.setAllDirichlet(); + else + values.setAllNeumann(); + + } + + /*! + * \brief Evaluates the boundary conditions for a Dirichlet + * boundary segment + * + * \param values Stores the Dirichlet values for the conservation equations in + * \f$ [ \textnormal{unit of primary variable} ] \f$ + * \param globalPos The global position + */ + void dirichletAtPos(PrimaryVariables &values, const GlobalPosition &globalPos) const + { + initialAtPos(values, globalPos); + } + + /*! + * \brief Evaluates the boundary conditions for a Neumann boundary segment. + * + * \param values Stores the Neumann values for the conservation equations in + * \f$ [ \textnormal{unit of conserved quantity} / (m^(dim-1) \cdot s )] \f$ + * \param globalPos The globalPosition of the boundary interface + */ + void neumannAtPos(PrimaryVariables &values, + const GlobalPosition &globalPos) const + { + // initialize values to zero, i.e. no-flow Neumann boundary conditions + values = 0; + values = 0; + + if (this->timeManager().time() + this->timeManager().timeStepSize() < injectionDuration_ + && globalPos[1] < 15 + eps_ && globalPos[1] > 7 - eps_ && globalPos[0] > 0.9*this->bBoxMax()[0]) + { + // inject air. negative values mean injection + values[Indices::contiNEqIdx] = -1e-4; // kg/(s*m^2) + values[Indices::contiWEqIdx] = 0.0; + } + } + + // \} + + /*! + * \name Volume terms + */ + // \{ + + /*! + * \brief Evaluates the initial values for a control volume + * + * \param values Stores the initial values for the conservation equations in + * \f$ [ \textnormal{unit of primary variables} ] \f$ + * \param globalPos The global position + */ + void initialAtPos(PrimaryVariables &values, const GlobalPosition &globalPos) const + { + + // get the water density at atmospheric conditions + const Scalar densityW = FluidSystem::H2O::liquidDensity(temperature(), 1.0e5); + + // assume an intially hydrostatic liquid pressure profile + // note: we subtract rho_w*g*h because g is defined negative + const Scalar pw = 1.0e5 - densityW*this->gravity()[dimWorld-1]*(aquiferDepth_ - globalPos[dimWorld-1]); + + values[Indices::pressureIdx] = pw; + values[saturationIdx] = 0.0; + } + // \} + +private: + static constexpr Scalar eps_ = 1e-6; + std::string name_; //! Problem name + Scalar aquiferDepth_; //! Depth of the aquifer in m + Scalar injectionDuration_; //! Duration of the injection in seconds +}; +} //end namespace + +#endif diff --git a/tutorial/ex1/injection2pspatialparams.hh b/tutorial/ex1/injection2pspatialparams.hh new file mode 100644 index 0000000000000000000000000000000000000000..f97ecab3f0abdf2fb7005f7d7d3980939fd3009f --- /dev/null +++ b/tutorial/ex1/injection2pspatialparams.hh @@ -0,0 +1,249 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/***************************************************************************** + * See the file COPYING for full copying permissions. * + * * + * This program is free software: you can redistribute it and/or modify * + * it under the terms of the GNU General Public License as published by * + * the Free Software Foundation, either version 2 of the License, or * + * (at your option) any later version. * + * * + * This program is distributed in the hope that it will be useful, * + * but WITHOUT ANY WARRANTY; without even the implied warranty of * + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * + * GNU General Public License for more details. * + * * + * You should have received a copy of the GNU General Public License * + * along with this program. If not, see . * + *****************************************************************************/ +/*! + * \file + * + * \brief Definition of the spatial parameters for the injection problem + * which uses the isothermal two-phase two-component + * fully implicit model. + */ + +#ifndef DUMUX_INJECTION_SPATIAL_PARAMS_HH +#define DUMUX_INJECTION_SPATIAL_PARAMS_HH + +#include +#include +#include + +#include +#include + +#include + +namespace Dumux +{ + +//forward declaration +template +class InjectionSpatialParams; + +namespace Properties +{ +// The spatial parameters TypeTag +NEW_TYPE_TAG(InjectionSpatialParams); + +// Set the spatial parameters +SET_TYPE_PROP(InjectionSpatialParams, SpatialParams, InjectionSpatialParams); + +// Set the material law parameterized by absolute saturations +SET_TYPE_PROP(InjectionSpatialParams, + MaterialLaw, + EffToAbsLaw >); +} + +/*! + * \ingroup TwoPTwoCModel + * \ingroup ImplicitTestProblems + * \brief Definition of the spatial parameters for the injection problem + * which uses the isothermal two-phase two-component + * fully implicit model. + */ +template +class InjectionSpatialParams : public ImplicitSpatialParams +{ + typedef ImplicitSpatialParams ParentType; + typedef typename GET_PROP_TYPE(TypeTag, Grid) Grid; + typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView; + typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar; + typedef typename Grid::ctype CoordScalar; + typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem; + + enum { + dim=GridView::dimension, + dimWorld=GridView::dimensionworld + }; + + typedef Dune::FieldVector GlobalPosition; + + typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry; + typedef typename GridView::template Codim<0>::Entity Element; + +public: + typedef typename GET_PROP_TYPE(TypeTag, MaterialLaw) MaterialLaw; + typedef typename MaterialLaw::Params MaterialLawParams; + + /*! + * \brief The constructor + * + * \param gridView The grid view + */ + InjectionSpatialParams(const GridView &gridView) + : ParentType(gridView) + { + aquiferHeightFromBottom_ = 30.0; + + // intrinsic permeabilities + aquitardK_ = GET_RUNTIME_PARAM(TypeTag, Scalar, SpatialParams.PermeabilityAquitard); + aquiferK_ = GET_RUNTIME_PARAM(TypeTag, Scalar, SpatialParams.PermeabilityAquifer); + + // porosities + aquitardPorosity_ = 0.2; + aquiferPorosity_ = 0.4; + + // residual saturations + aquitardMaterialParams_.setSwr(0.2); + aquitardMaterialParams_.setSnr(0.0); + aquiferMaterialParams_.setSwr(0.2); + aquiferMaterialParams_.setSnr(0.0); + + // parameters for the Brooks-Corey law + aquitardMaterialParams_.setPe(GET_RUNTIME_PARAM(TypeTag, Scalar, SpatialParams.EntryPressureAquitard)); + aquiferMaterialParams_.setPe(GET_RUNTIME_PARAM(TypeTag, Scalar, SpatialParams.EntryPressureAquifer)); + aquitardMaterialParams_.setLambda(2.0); + aquiferMaterialParams_.setLambda(2.0); + + } + + /*! + * \brief Returns the intrinsic permeability tensor \f$[m^2]\f$ + * + * \param element The finite element + * \param fvGeometry The finite volume geometry of the element + * \param scvIdx The local index of the sub-control volume + */ + const Scalar intrinsicPermeability(const Element &element, + const FVElementGeometry &fvGeometry, + const int scvIdx) const + { + const GlobalPosition &globalPos = fvGeometry.subContVol[scvIdx].global; + if (isInAquitard_(globalPos)) + return aquitardK_; + return aquiferK_; + } + + /*! + * \brief Returns the porosity \f$[-]\f$ + * + * \param element The finite element + * \param fvGeometry The finite volume geometry of the element + * \param scvIdx The local index of the sub-control volume + */ + Scalar porosity(const Element &element, + const FVElementGeometry &fvGeometry, + const int scvIdx) const + { + const GlobalPosition &globalPos = fvGeometry.subContVol[scvIdx].global; + if (isInAquitard_(globalPos)) + return aquitardPorosity_; + return aquiferPorosity_; + } + + + /*! + * \brief Returns the parameter object for the capillary-pressure/ + * saturation material law + * + * \param element The finite element + * \param fvGeometry The finite volume geometry of the element + * \param scvIdx The local index of the sub-control volume + */ + const MaterialLawParams& materialLawParams(const Element &element, + const FVElementGeometry &fvGeometry, + const int scvIdx) const + { + const GlobalPosition &globalPos = fvGeometry.subContVol[scvIdx].global; + if (isInAquitard_(globalPos)) + return aquitardMaterialParams_; + return aquiferMaterialParams_; + } + + /*! + * These parameters are only needed for nonisothermal models. Comment them in if you want to implement the 2pni model. + */ + + /*! + * \brief Returns the heat capacity \f$[J / (kg K)]\f$ of the rock matrix. + * + * This is only required for non-isothermal models. + * + * \param element The finite element + * \param fvGeometry The finite volume geometry + * \param scvIdx The local index of the sub-control volume + */ +// Scalar solidHeatCapacity(const Element &element, +// const FVElementGeometry &fvGeometry, +// const int scvIdx) const +// { +// return 790; // specific heat capacity of granite [J / (kg K)] +// } + + /*! + * \brief Returns the mass density \f$[kg / m^3]\f$ of the rock matrix. + * + * This is only required for non-isothermal models. + * + * \param element The finite element + * \param fvGeometry The finite volume geometry + * \param scvIdx The local index of the sub-control volume + */ +// Scalar solidDensity(const Element &element, +// const FVElementGeometry &fvGeometry, +// const int scvIdx) const +// { +// return 2700; // density of granite [kg/m^3] +// } + + /*! + * \brief Returns the thermal conductivity \f$\mathrm{[W/(m K)]}\f$ of the solid + * + * This is only required for non-isothermal models. + * + * \param element The finite element + * \param fvGeometry The finite volume geometry of the element + * \param scvIdx The local index of the sub-control volume + */ +// Scalar solidThermalConductivity(const Element &element, +// const FVElementGeometry &fvGeometry, +// const int scvIdx) const +// { +// return lambdaSolid_; +// } + +private: + + static constexpr Scalar eps_ = 1e-6; + + bool isInAquitard_(const GlobalPosition &globalPos) const + { return globalPos[dimWorld-1] > aquiferHeightFromBottom_ + eps_; } + + Scalar aquitardK_; + Scalar aquiferK_; + Scalar aquiferHeightFromBottom_; + + + Scalar aquitardPorosity_; + Scalar aquiferPorosity_; + + MaterialLawParams aquitardMaterialParams_; + MaterialLawParams aquiferMaterialParams_; +}; + +} + +#endif diff --git a/tutorial/ex2/CMakeLists.txt b/tutorial/ex2/CMakeLists.txt new file mode 100644 index 0000000000000000000000000000000000000000..d53f775c2a43deb9104ec172710926fa86af0a7e --- /dev/null +++ b/tutorial/ex2/CMakeLists.txt @@ -0,0 +1,6 @@ +# a compositional two-phase simulation program +dune_add_test(NAME exercise2 + SOURCES exercise2.cc) + +# add a symlink for the input file +dune_symlink_to_source_files(FILES "exercise2.input") diff --git a/tutorial/ex2/README.md b/tutorial/ex2/README.md new file mode 100644 index 0000000000000000000000000000000000000000..5473aaba6bb8d65865568721865b254e12e1a348 --- /dev/null +++ b/tutorial/ex2/README.md @@ -0,0 +1,200 @@ +# Exercise #2 (DuMuX course) +
+## Problem set-up + +The problem setup is identical to the previous [_exercise 1_](../ex1/README.md). + +## Preparing the exercise + +* Navigate to the directory `dumux/tutorial/ex2` + +_Exercise 2_ deals with a two-phase compositional problem (__2p2c__). Goal is to learn how to use compile and runtime parameters and the _DuMuX property system_. + +

+### Task 1: Getting familiar with the code +
+ +Locate all the files you will need for this exercise +* The __main file__: `exercise2.cc` +* The __problem file__: `injection2p2cproblem.hh` +* The __spatial parameters file__: `injection2p2cspatialparams.hh` +* The __input file__: `exercise2.input` +* Two header files containing: + * a custom __local residual__ in: `mylocalresidual.hh` + * a custom __material law__ in: `mymateriallaw.hh` + + +


+### Task 2: Compiling and running the program +
+ +* Change to the build-directory + +```bash +cd ../../build-cmake/tutorial/exercise2 +``` + +* Compile the executable `exercise2` + +```bash +make exercise2 +``` + +* Run the problem and inspect the result + +```bash +./exercise2 +``` +Note: Because the input file has the same name as the +executable, DuMuX will find it automatically. + +If gnuplot is installed on your system, you should see a plot of the capillary pressure - saturation relationship. + +


+### Task 3: Implement and use a different material law +
+ +DuMuX uses the term _material law_ to describe the law used to compute +* pc-Sw relations +* kr-Sw relations +* their inverse relations + +The file `mymateriallaw.hh` contains a custom implementation of such a material law. + +* Implement the method `Scalar pc(Scalar sw)` by implementing your own capillary pressure relationship, e.g. a simple linear relationship $`p_C(S_w) = 1\cdot 10^5 \cdot (1-S_w) + p_e`$. + +Note: `MyMaterialLaw` uses the `BrooksCoreyParams` class as parameter input. You can get the entry pressure that is set in the spatial params as follows + +```c++ +const auto pe = params.pe(); +``` + +The type (i.e. C++ type) of the material law is set in the file `injection2pspatialparams.hh` by using the DuMuX property system + +```c++ +SET_PROP(InjectionSpatialParams, MaterialLaw) +{ + using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar); + using type = EffToAbsLaw>; +}; +``` + +* Make DuMuX use your own material law by including the header `mymateriallaw.hh` and changing the alias `type`. This will make sure that your material law is used everywhere else in the code. + +Note: Also use the wrapper class `EffToAbsLaw`. It takes care of converting absolute to effective saturations considering residual saturations. `MyMaterialLaw` +as other material laws (like Brooks-Corey, VanGenuchten, ...) in DuMuX only deals with effective saturations. + +* Verify your changes by recompiling and running the program. You should see a plot of your new function. + +For the next task, disable the plotting feature by changing the settings in the input file `exercise2.input` + +```ini +[Problem] +OnlyPlotMaterialLaws = false +``` + +


+### Task 4: Enable/Disable Gravity -> DuMuX parameters +
+ +DuMuX has many parameters that have default values. For example, all simulations consider gravity effects by default. +You can disable gravity for a study, simply by setting the parameter in the input file + +```ini +[Problem] +EnableGravity = false +``` + +Run the simulation with and without gravity. Change the `Problem.Name` parameter to create output files with different +names. Compare the results using paraview. You should immediately see the difference. + +A list of parameters that can be set through the input file is given [here](http://www.dumux.org/doxygen-stable/html-2.11/a06387.php). + + +


+### Task 5: Implement your own local residual +
+ +Most types in DuMuX are properties that can be changed just like the material law. In the following task we implement our own 2p2c local residual, i.e. the class that computes the element residual in every Newton step. The file `mylocalresidual.hh` contains a copy of the original local residual class used for the 2p2c model renamed to `template class MyTwoPTwoCLocalResidual`. + +* Make DuMuX use this new local residual by inluding the header `mylocalresidual.hh` and setting the corresponding property in the `Property` namespace in the file `injection2p2cproblem.hh` + +```c++ +// note that every property struct knows about TypeTag +SET_PROP(Injection2p2cProblem, LocalResidual) +{ + using type = MyTwoPTwoCLocalResidualocal; +}; + +// or using the convenience macro +SET_TYPE_PROP(Injection2p2cProblem, LocalResidual, + MyTwoPTwoCLocalResidualocal); +``` + +* Implement an output to the terminal in the constructor of `MyTwoPTwoCLocalResidual` e.g. + +```c++ +MyTwoPTwoCLocalResidual() +{ + std::cout << "Using MyTwoPTwoCLocalResidual." << std::endl; +} +``` + +* Verify you are using the new class by compiling and running the new program and inspecting the terminal output. + +You want to make the new local residual special by adding a switch enabling / disabling diffusion. We will achieve this with a DuMuX parameter which is read from the input file and defaults to a property value if the input file doesn't contain the parameter. + +* Create a new `TypeTag` node, a new `PropertyTag`, and set a default in the `mylocalresidual.hh` file by adding + +```c++ +namespace Dumux { + +namespace Properties +{ + NEW_TYPE_TAG(MyLocalResidualParams); // creates a new TypeTag node + NEW_PROP_TAG(ProblemEnableDiffusion); // creates a new property + SET_BOOL_PROP(MyLocalResidualParams, + ProblemEnableDiffusion, true); // set a default value +} +... +``` + +* Then enhance the problem TypeTag (the master TypeTag) by inheriting from the created node (in `injection2p2cproblem.hh`) + +```c++ +... + NEW_TYPE_TAG(Injection2p2cProblem, // the problem TypeTag + INHERITS_FROM(TwoPTwoC, // the 2p2c model TypeTag node + InjectionSpatialParams, // the spatial params TypeTag node + MyLocalResidualParams)); // the local residual params TypeTag node +... +``` +Note: the property inheritance graph now looks like this + +![property tree](../extradoc/exercise2_properties.png) + + +* Modify the `computeFlux` method to only call the `diffusiveFlux` method if diffusion is enabled. You can get the new parameter by adding the lines + +```c++ +// ... in the constructor of MyTwoPTwoCLocalResidual + enableDiffusion_ = GET_PARAM_FROM_GROUP(TypeTag, bool, + Problem, EnableDiffusion); + +// ... in the private member section of MyTwoPTwoCLocalResidual +private: + bool enableDiffusion_; +``` + +You can now enable and disable diffusion through the input file + +```ini +[Problem] +EnableDiffusion = true / false +``` + +* Verify the difference in the parameter $`x_w^{N2}`$, i.e. the mole fraction of nitrogen in the +water phase, with and without diffusion. + +Note that due to diffusion being a slow process you +can only see the difference in later times. diff --git a/tutorial/ex2/exercise2.cc b/tutorial/ex2/exercise2.cc new file mode 100644 index 0000000000000000000000000000000000000000..aba0b0e8c76986dd481d8ad9978ddb7874cb01ae --- /dev/null +++ b/tutorial/ex2/exercise2.cc @@ -0,0 +1,67 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/***************************************************************************** + * See the file COPYING for full copying permissions. * + * * + * This program is free software: you can redistribute it and/or modify * + * it under the terms of the GNU General Public License as published by * + * the Free Software Foundation, either version 2 of the License, or * + * (at your option) any later version. * + * * + * This program is distributed in the hope that it will be useful, * + * but WITHOUT ANY WARRANTY; without even the implied warranty of * + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * + * GNU General Public License for more details. * + * * + * You should have received a copy of the GNU General Public License * + * along with this program. If not, see . * + *****************************************************************************/ +/*! + * \file + * + * \brief Test for the two-phase two-component CC model. + */ +#include +#include "injection2p2cproblem.hh" +#include + +/*! + * \brief Provides an interface for customizing error messages associated with + * reading in parameters. + * + * \param progName The name of the program, that was tried to be started. + * \param errorMsg The error message that was issued by the start function. + * Comprises the thing that went wrong and a general help message. + */ +void usage(const char *progName, const std::string &errorMsg) +{ + if (errorMsg.size() > 0) { + std::string errorMessageOut = "\nUsage: "; + errorMessageOut += progName; + errorMessageOut += " [options]\n"; + errorMessageOut += errorMsg; + errorMessageOut += "\n\nThe list of mandatory options for this program is:\n" + "\t-TimeManager.TEnd End of the simulation [s] \n" + "\t-TimeManager.DtInitial Initial timestep size [s] \n" + "\t-Grid.File Name of the file containing the grid \n" + "\t definition in DGF format\n" + "\t-FluidSystem.NTemperature Number of tabularization entries [-] \n" + "\t-FluidSystem.NPressure Number of tabularization entries [-] \n" + "\t-FluidSystem.PressureLow Low end for tabularization of fluid properties [Pa] \n" + "\t-FluidSystem.PressureHigh High end for tabularization of fluid properties [Pa] \n" + "\t-FluidSystem.TemperatureLow Low end for tabularization of fluid properties [Pa] \n" + "\t-FluidSystem.TemperatureHigh High end for tabularization of fluid properties [Pa] \n" + "\t-SimulationControl.Name The name of the output files [-] \n" + "\t-InitialConditions.Temperature Initial temperature in the reservoir [K] \n" + "\t-InitialConditions.DepthBOR Depth below ground surface [m] \n"; + + std::cout << errorMessageOut + << "\n"; + } +} + +int main(int argc, char** argv) +{ + typedef TTAG(Injection2p2pcCCProblem) ProblemTypeTag; + return Dumux::start(argc, argv, usage); +} diff --git a/tutorial/ex2/exercise2.input b/tutorial/ex2/exercise2.input new file mode 100644 index 0000000000000000000000000000000000000000..c7149509ad8f494364f74f05a70e66b935f51c52 --- /dev/null +++ b/tutorial/ex2/exercise2.input @@ -0,0 +1,25 @@ +[TimeManager] +DtInitial = 3600 # in seconds +TEnd = 3.154e9 # in seconds, i.e ten years + +[Grid] +LowerLeft = 0 0 +UpperRight = 60 40 +Cells = 24 16 + +[Problem] +Name = infiltration +OnlyPlotMaterialLaws = true +AquiferDepth = 2700.0 # m +TotalAreaSpecificInflow = 1e-4 # kg / (s*m^2) +InjectionDuration = 2.628e6 # in seconds, i.e. one month + +# TODO: dumux-course-task +# Set Problem.EnableGravity +# Set Problem.EnableDiffusion + +[SpatialParams] +PermeabilityAquitard = 1e-15 # m^2 +EntryPressureAquitard = 4.5e4 # Pa +PermeabilityAquifer = 1e-12 # m^2 +EntryPressureAquifer = 1e4 # Pa diff --git a/tutorial/ex2/injection2p2cproblem.hh b/tutorial/ex2/injection2p2cproblem.hh new file mode 100644 index 0000000000000000000000000000000000000000..6b611593e14ed25de875bcf24d9ec1efe3e92737 --- /dev/null +++ b/tutorial/ex2/injection2p2cproblem.hh @@ -0,0 +1,316 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/***************************************************************************** + * See the file COPYING for full copying permissions. * + * * + * This program is free software: you can redistribute it and/or modify * + * it under the terms of the GNU General Public License as published by * + * the Free Software Foundation, either version 2 of the License, or * + * (at your option) any later version. * + * * + * This program is distributed in the hope that it will be useful, * + * but WITHOUT ANY WARRANTY; without even the implied warranty of * + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * + * GNU General Public License for more details. * + * * + * You should have received a copy of the GNU General Public License * + * along with this program. If not, see . * + *****************************************************************************/ +/*! + * \file + * + * \brief Problem where air is injected under a low permeable layer in a depth of 2700m. + */ +#ifndef DUMUX_INJECTION_2P2C_PROBLEM_HH +#define DUMUX_INJECTION_2P2C_PROBLEM_HH + +#include +#include +#include + +#include "injection2p2cspatialparams.hh" + +// TODO: dumux-course-task +// Include the local residual header + +namespace Dumux +{ + +// foward declaration +template +class Injection2p2cProblem; + +// setup property TypeTag +namespace Properties +{ +// TODO: dumux-course-task +// inherit from MyLocalResidualParams +NEW_TYPE_TAG(Injection2p2cProblem, INHERITS_FROM(TwoPTwoC, InjectionSpatialParams)); +NEW_TYPE_TAG(Injection2p2cBoxProblem, INHERITS_FROM(BoxModel, Injection2p2cProblem)); +NEW_TYPE_TAG(Injection2p2pcCCProblem, INHERITS_FROM(CCModel, Injection2p2cProblem)); + +// Set the grid type +SET_TYPE_PROP(Injection2p2cProblem, Grid, Dune::YaspGrid<2>); + +// Set the problem property +SET_TYPE_PROP(Injection2p2cProblem, Problem, Injection2p2cProblem); + +// TODO: dumux-course-task +// change the local residual type to MyTwoPTwoCLocalResidual + +// Set fluid configuration +SET_TYPE_PROP(Injection2p2cProblem, + FluidSystem, + FluidSystems::H2ON2); + +// Define whether mole(true) or mass (false) fractions are used +SET_BOOL_PROP(Injection2p2cProblem, UseMoles, true); + +} // end namespace Properties + +/*! + * \ingroup TwoPTwoCModel + * \ingroup ImplicitTestProblems + * \brief Problem where air is injected under a low permeable layer in a depth of 2700m. + * + * The domain is sized 60m times 40m and consists of two layers, a moderately + * permeable one for \f$ y<22m\f$ and one with a lower permeablility + * in the rest of the domain. + * + * Nitrogen is injected into a water-filled aquifer through a well. First, we inject for one month. + * Then, we continue simulating the development of the nitrogen plume for 10 years. + * This is realized with a Neumann boundary condition at the right boundary + * (\f$ 7m +class Injection2p2cProblem : public ImplicitPorousMediaProblem +{ + typedef ImplicitPorousMediaProblem ParentType; + + typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView; + typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar; + typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem; + + // grid world dimension + static constexpr auto dimWorld = GridView::dimensionworld; + + typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices; + typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables; + typedef typename GET_PROP_TYPE(TypeTag, BoundaryTypes) BoundaryTypes; + typedef typename GET_PROP_TYPE(TypeTag, TimeManager) TimeManager; + typedef typename GET_PROP_TYPE(TypeTag, MaterialLaw) MaterialLaw; + typedef Dune::FieldVector GlobalPosition; + +public: + /*! + * \brief The constructor + * + * \param timeManager The time manager + * \param gridView The grid view + */ + Injection2p2cProblem(TimeManager &timeManager, const GridView &gridView) + : ParentType(timeManager, gridView) + { + // initialize the tables of the fluid system + FluidSystem::init(/*tempMin=*/273.15, + /*tempMax=*/423.15, + /*numTemp=*/50, + /*pMin=*/0.0, + /*pMax=*/30e6, + /*numP=*/300); + + // name of the problem and output file + name_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, std::string, Problem, Name); + // depth of the aquifer, units: m + aquiferDepth_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Problem, AquiferDepth); + // inflow rate of nitrogen water vapor mixture, units: kg/(s m^2) + totalAreaSpecificInflow_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Problem, TotalAreaSpecificInflow); + // the duration of the injection, units: second + injectionDuration_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Problem, InjectionDuration); + } + + /*! + * \brief User defined output after the time integration + * + * Will be called diretly after the time integration + */ + void postTimeStep() + { + // Calculate storage terms + PrimaryVariables storageW, storageN; + this->model().globalPhaseStorage(storageW, Indices::wPhaseIdx); + this->model().globalPhaseStorage(storageN, Indices::nPhaseIdx); + + // Write mass balance information for rank 0 + if (this->gridView().comm().rank() == 0) + { + std::cout <<"Storage: wetting=[" << storageW << "]" + << " nonwetting=[" << storageN << "]" << std::endl; + } + } + + /*! + * \name Problem parameters + */ + // \{ + + /*! + * \brief Returns the problem name + * + * This is used as a prefix for files generated by the simulation. + */ + const std::string& name() const + { return name_; } + + /*! + * \brief Returns the temperature in \f$ K \f$ + */ + Scalar temperature() const + { return 273.15 + 30; } + + /*! + * \brief Returns the source term + * + * \param values Stores the source values for the conservation equations in + * \f$ [ \textnormal{unit of primary variable} / (m^\textrm{dim} \cdot s )] \f$ + * \param globalPos The global position + */ + void sourceAtPos(PrimaryVariables &values, + const GlobalPosition &globalPos) const + { values = 0; } + + // \} + + /*! + * \name Boundary conditions + */ + // \{ + + /*! + * \brief Specifies which kind of boundary condition should be + * used for which equation on a given boundary segment + * + * \param values Stores the value of the boundary type + * \param globalPos The global position + */ + void boundaryTypesAtPos(BoundaryTypes &values, + const GlobalPosition &globalPos) const + { + // Set Dirichlet at the bottom of the domain + if (globalPos[dimWorld-1] < eps_) + { + values.setAllDirichlet(); + } + + // and Neuman boundary conditions everywhere else + // note that we don't differentiate between Neumann and Robin boundary types + else + { + values.setAllNeumann(); + } + } + + /*! + * \brief Evaluates the boundary conditions for a Dirichlet + * boundary segment + * + * \param values Stores the Dirichlet values for the conservation equations in + * \f$ [ \textnormal{unit of primary variable} ] \f$ + * \param globalPos The global position + */ + void dirichletAtPos(PrimaryVariables &values, const GlobalPosition &globalPos) const + { initialAtPos(values, globalPos); } + + /*! + * \brief Evaluates the boundary conditions for a Neumann boundary segment. + * + * \param values Stores the Neumann values for the conservation equations in + * \f$ [ \textnormal{unit of conserved quantity} / (m^(dim-1) \cdot s )] \f$ + * \param globalPos The globalPosition of the boundary interface + */ + void neumannAtPos(PrimaryVariables &values, + const GlobalPosition &globalPos) const + { + // initialize values to zero, i.e. no-flow Neumann boundary conditions + values = 0; + + //if we are inside the injection zone set inflow Neumann boundary conditions + if (this->timeManager().time() + this->timeManager().timeStepSize() < injectionDuration_ + && globalPos[1] < 15 + eps_ && globalPos[1] > 7 - eps_ && globalPos[0] > 0.9*this->bBoxMax()[0]) + { + // set the Neumann values for the Nitrogen component balance + // convert from units kg/(s*m^2) to mole/(s*m^2) + values[Indices::contiNEqIdx] = -totalAreaSpecificInflow_/FluidSystem::molarMass(FluidSystem::nCompIdx); + values[Indices::contiWEqIdx] = 0.0; + } + } + + // \} + + /*! + * \name Volume terms + */ + // \{ + + /*! + * \brief Evaluates the initial values for a control volume + * + * \param values Stores the initial values for the conservation equations in + * \f$ [ \textnormal{unit of primary variables} ] \f$ + * \param globalPos The global position + */ + void initialAtPos(PrimaryVariables &values, const GlobalPosition &globalPos) const + { + // get the water density at atmospheric conditions + const Scalar densityW = FluidSystem::H2O::liquidDensity(temperature(), 1.0e5); + + // assume an intially hydrostatic liquid pressure profile + // note: we subtract rho_w*g*h because g is defined negative + const Scalar pw = 1.0e5 - densityW*this->gravity()[dimWorld-1]*(aquiferDepth_ - globalPos[dimWorld-1]); + + // initially we have some nitrogen dissolved + // saturation mole fraction would be + // moleFracLiquidN2 = (pw + pc + p_vap^sat)/henry; + const Scalar moleFracLiquidN2 = pw*0.95/BinaryCoeff::H2O_N2::henry(temperature()); + + // note that because we start with a single phase system the primary variables + // are pl and x^w_N2. This will switch as soon after we start injecting to a two + // phase system so the primary variables will be pl and Sn (non-wetting saturation). + values[Indices::switchIdx] = moleFracLiquidN2; + values[Indices::pressureIdx] = pw; + } + + /*! + * \brief Return the initial phase state inside a control volume. + * + * \param globalPos The global position + * \note we start with a single phase system + */ + int initialPhasePresenceAtPos(const GlobalPosition &globalPos) const + { return Indices::wPhaseOnly; } + + // \} + + //! If we should write restart files + bool shouldWriteRestartFile() const + { return false; } + +private: + static constexpr Scalar eps_ = 1e-6; + std::string name_; //! Problem name + Scalar aquiferDepth_; //! Depth of the aquifer in m + Scalar totalAreaSpecificInflow_; //! Area specific inflow rate in mole/(s*m^2) + Scalar injectionDuration_; //! Duration of the injection in seconds +}; + +} // end namespace Dumux + +#endif diff --git a/tutorial/ex2/injection2p2cspatialparams.hh b/tutorial/ex2/injection2p2cspatialparams.hh new file mode 100644 index 0000000000000000000000000000000000000000..aac0effb442c7ecf28ad030c74bc4e18db14c523 --- /dev/null +++ b/tutorial/ex2/injection2p2cspatialparams.hh @@ -0,0 +1,221 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/***************************************************************************** + * See the file COPYING for full copying permissions. * + * * + * This program is free software: you can redistribute it and/or modify * + * it under the terms of the GNU General Public License as published by * + * the Free Software Foundation, either version 2 of the License, or * + * (at your option) any later version. * + * * + * This program is distributed in the hope that it will be useful, * + * but WITHOUT ANY WARRANTY; without even the implied warranty of * + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * + * GNU General Public License for more details. * + * * + * You should have received a copy of the GNU General Public License * + * along with this program. If not, see . * + *****************************************************************************/ +/*! + * \file + * + * \brief Definition of the spatial parameters for the injection problem + * which uses the isothermal two-phase two-component + * fully implicit model. + */ + +#ifndef DUMUX_INJECTION_SPATIAL_PARAMS_HH +#define DUMUX_INJECTION_SPATIAL_PARAMS_HH + +#include +#include +#include +// TODO: dumux-course-task +// Inlcude your own material law + +#include +#include + +#include + +namespace Dumux +{ + +// forward declaration +template +class InjectionSpatialParams; + +// setup property TypeTag +namespace Properties +{ +// The spatial parameters TypeTag +NEW_TYPE_TAG(InjectionSpatialParams); + +// Set the spatial parameters +SET_TYPE_PROP(InjectionSpatialParams, SpatialParams, InjectionSpatialParams); + +// TODO: dumux-course-task +// Use your own material law instead +// Set the material law parameterized by absolute saturations +SET_PROP(InjectionSpatialParams, MaterialLaw) +{ + using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar); + using type = EffToAbsLaw>; +}; + +} // end namespace Properties + +/*! + * \ingroup TwoPTwoCModel + * \ingroup ImplicitTestProblems + * \brief Definition of the spatial parameters for the injection problem + * which uses the isothermal two-phase two-component + * fully implicit model. + */ +template +class InjectionSpatialParams : public ImplicitSpatialParams +{ + typedef ImplicitSpatialParams ParentType; + typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView; + typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar; + typedef typename GridView::ctype CoordinateType; + + enum { + dim=GridView::dimension, + dimWorld=GridView::dimensionworld + }; + + typedef Dune::FieldVector GlobalPosition; + typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry; + typedef typename GridView::template Codim<0>::Entity Element; + typedef typename GET_PROP_TYPE(TypeTag, MaterialLaw) MaterialLaw; + typedef typename MaterialLaw::Params MaterialLawParams; + +public: + + /*! + * \brief The constructor + * + * \param gridView The grid view + */ + InjectionSpatialParams(const GridView &gridView) + : ParentType(gridView) + { + aquiferHeightFromBottom_ = 30.0; + + // intrinsic permeabilities + aquitardK_ = GET_RUNTIME_PARAM(TypeTag, Scalar, SpatialParams.PermeabilityAquitard); + aquiferK_ = GET_RUNTIME_PARAM(TypeTag, Scalar, SpatialParams.PermeabilityAquifer); + + // porosities + aquitardPorosity_ = 0.2; + aquiferPorosity_ = 0.4; + + // residual saturations + aquitardMaterialParams_.setSwr(0.2); + aquitardMaterialParams_.setSnr(0.0); + aquiferMaterialParams_.setSwr(0.2); + aquiferMaterialParams_.setSnr(0.0); + + // parameters for the Brooks-Corey law + aquitardMaterialParams_.setPe(GET_RUNTIME_PARAM(TypeTag, Scalar, SpatialParams.EntryPressureAquitard)); + aquiferMaterialParams_.setPe(GET_RUNTIME_PARAM(TypeTag, Scalar, SpatialParams.EntryPressureAquifer)); + aquitardMaterialParams_.setLambda(2.0); + aquiferMaterialParams_.setLambda(2.0); + + // plot the material laws using gnuplot and exit + if (GET_RUNTIME_PARAM(TypeTag, bool, Problem.OnlyPlotMaterialLaws)) + { + plotMaterialLaws(); + exit(0); + } + } + + /*! + * \brief Returns the intrinsic permeability tensor \f$[m^2]\f$ + * + * \param element The finite element + * \param fvGeometry The finite volume geometry of the element + * \param scvIdx The local index of the sub-control volume + */ + Scalar intrinsicPermeability(const Element &element, + const FVElementGeometry &fvGeometry, + const int scvIdx) const + { + const GlobalPosition &globalPos = fvGeometry.subContVol[scvIdx].global; + if (isInAquitard_(globalPos)) + return aquitardK_; + return aquiferK_; + } + + /*! + * \brief Returns the porosity \f$[-]\f$ + * + * \param element The finite element + * \param fvGeometry The finite volume geometry of the element + * \param scvIdx The local index of the sub-control volume + */ + Scalar porosity(const Element &element, + const FVElementGeometry &fvGeometry, + const int scvIdx) const + { + const GlobalPosition &globalPos = fvGeometry.subContVol[scvIdx].global; + if (isInAquitard_(globalPos)) + return aquitardPorosity_; + return aquiferPorosity_; + } + + + /*! + * \brief Returns the parameter object for the capillary-pressure/ + * saturation material law + * + * \param element The finite element + * \param fvGeometry The finite volume geometry of the element + * \param scvIdx The local index of the sub-control volume + */ + const MaterialLawParams& materialLawParams(const Element &element, + const FVElementGeometry &fvGeometry, + const int scvIdx) const + { + const GlobalPosition &globalPos = fvGeometry.subContVol[scvIdx].global; + if (isInAquitard_(globalPos)) + return aquitardMaterialParams_; + return aquiferMaterialParams_; + } + + /*! + * \brief Creates a gnuplot output of the pc-Sw curve + */ + void plotMaterialLaws() + { + PlotMaterialLaw plotMaterialLaw; + GnuplotInterface gnuplot; + plotMaterialLaw.addpcswcurve(gnuplot, aquitardMaterialParams_, 0.2, 1.0, "upper layer (fine, aquitard)", "w lp"); + plotMaterialLaw.addpcswcurve(gnuplot, aquiferMaterialParams_, 0.2, 1.0, "lower layer (coarse, aquifer)", "w l"); + gnuplot.setOption("set xrange [0:1]"); + gnuplot.setOption("set label \"residual\\nsaturation\" at 0.1,100000 center"); + gnuplot.plot("pc-Sw"); + } + +private: + + static constexpr Scalar eps_ = 1e-6; + + bool isInAquitard_(const GlobalPosition &globalPos) const + { return globalPos[dimWorld-1] > aquiferHeightFromBottom_ + eps_; } + + Scalar aquitardK_; + Scalar aquiferK_; + Scalar aquiferHeightFromBottom_; + + Scalar aquitardPorosity_; + Scalar aquiferPorosity_; + + MaterialLawParams aquitardMaterialParams_; + MaterialLawParams aquiferMaterialParams_; +}; + +} // end namespace Dumux + +#endif diff --git a/tutorial/ex2/mylocalresidual.hh b/tutorial/ex2/mylocalresidual.hh new file mode 100644 index 0000000000000000000000000000000000000000..bcaae17364c68a0498f8e1cd6f16232883f63649 --- /dev/null +++ b/tutorial/ex2/mylocalresidual.hh @@ -0,0 +1,535 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/***************************************************************************** + * See the file COPYING for full copying permissions. * + * * + * This program is free software: you can redistribute it and/or modify * + * it under the terms of the GNU General Public License as published by * + * the Free Software Foundation, either version 2 of the License, or * + * (at your option) any later version. * + * * + * This program is distributed in the hope that it will be useful, * + * but WITHOUT ANY WARRANTY; without even the implied warranty of * + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * + * GNU General Public License for more details. * + * * + * You should have received a copy of the GNU General Public License * + * along with this program. If not, see . * + *****************************************************************************/ +/*! + * \file + * + * \brief Element-wise calculation of the Jacobian matrix for problems + * using the two-phase two-component fully implicit model. + */ + +#ifndef DUMUX_MY_2P2C_LOCAL_RESIDUAL_HH +#define DUMUX_MY_2P2C_LOCAL_RESIDUAL_HH + +#include + +namespace Dumux +{ + +// TODO: dumux-course-task +// implement new TypeTag MyLocalResidualParams + +/*! + * \ingroup TwoPTwoCModel + * \ingroup ImplicitLocalResidual + * \brief Element-wise calculation of the Jacobian matrix for problems + * using the two-phase two-component fully implicit model. + * + * This class is used to fill the gaps in ImplicitLocalResidual for the + * two-phase two-component flow. + */ +template +class MyTwoPTwoCLocalResidual: public GET_PROP_TYPE(TypeTag, BaseLocalResidual) +{ + protected: + typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar; + typedef typename GET_PROP_TYPE(TypeTag, LocalResidual) Implementation; + typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem; + typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables; + typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry; + typedef typename GET_PROP_TYPE(TypeTag, VolumeVariables) VolumeVariables; + typedef typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables) ElementVolumeVariables; + typedef typename GET_PROP_TYPE(TypeTag, ElementBoundaryTypes) ElementBoundaryTypes; + typedef typename GET_PROP_TYPE(TypeTag, FluxVariables) FluxVariables; + enum + { + numPhases = GET_PROP_VALUE(TypeTag, NumPhases), + numComponents = GET_PROP_VALUE(TypeTag, NumComponents) + }; + + typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices; + enum + { + contiWEqIdx = Indices::contiWEqIdx, + contiNEqIdx = Indices::contiNEqIdx, + wPhaseIdx = Indices::wPhaseIdx, + nPhaseIdx = Indices::nPhaseIdx, + wCompIdx = Indices::wCompIdx, + nCompIdx = Indices::nCompIdx + }; + + typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView; + typedef typename GridView::template Codim<0>::Entity Element; + + static constexpr unsigned int replaceCompEqIdx = GET_PROP_VALUE(TypeTag, ReplaceCompEqIdx); + + //! Property that defines whether mole or mass fractions are used + static const bool useMoles = GET_PROP_VALUE(TypeTag, UseMoles); + + public: + /*! + * \brief Constructor + * + * Sets the mass upwind weight. + */ + MyTwoPTwoCLocalResidual() + { + // retrieve the upwind weight for the mass conservation equations. Use the value + // specified via the property system as default, and overwrite + // it by the run-time parameter from the Dune::ParameterTree + massUpwindWeight_ = GET_PARAM_FROM_GROUP(TypeTag, Scalar, Implicit, MassUpwindWeight); + + // TODO: dumux-course-task + // get parameter Problem.EnableDiffusion + } + + /*! + * \brief Evaluate the storage term of the current solution in a + * single phase. + * + * \param element The element + * \param phaseIdx The index of the fluid phase + */ + void evalPhaseStorage(const Element &element, const int phaseIdx) + { + FVElementGeometry fvGeometry; + fvGeometry.update(this->gridView_(), element); + ElementBoundaryTypes bcTypes; + bcTypes.update(this->problem_(), element, fvGeometry); + ElementVolumeVariables elemVolVars; + elemVolVars.update(this->problem_(), element, fvGeometry, false); + + this->storageTerm_.resize(fvGeometry.numScv); + this->storageTerm_ = 0; + + this->elemPtr_ = &element; + this->fvElemGeomPtr_ = &fvGeometry; + this->bcTypesPtr_ = &bcTypes; + this->prevVolVarsPtr_ = 0; + this->curVolVarsPtr_ = &elemVolVars; + evalPhaseStorage_(phaseIdx); + } + + /*! + * \brief Evaluate the amount of all conservation quantities + * (e.g. phase mass) within a sub-control volume. + * + * \param storage The mass of the component within the sub-control volume + * \param scvIdx The sub-control-volume index + * \param usePrevSol Based on usePrevSol solution of current or previous time step is used + * + * The result should be averaged over the volume (e.g. phase mass + * inside a sub-control volume divided by the volume) + */ + void computeStorage(PrimaryVariables &storage, const int scvIdx, bool usePrevSol) const + { + // if flag usePrevSol is set, the solution from the previous + // time step is used, otherwise the current solution is + // used. The secondary variables are used accordingly. This + // is required to compute the derivative of the storage term + // using the implicit Euler method. + const ElementVolumeVariables &elemVolVars = usePrevSol ? this->prevVolVars_() + : this->curVolVars_(); + const VolumeVariables &volVars = elemVolVars[scvIdx]; + + // compute storage term of all components within all phases + storage = 0; + if(useMoles) // mole-fraction formulation + { + for (unsigned int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) + { + for (unsigned int compIdx = contiCompIdx1_(); compIdx <= contiCompIdx2_(); ++compIdx) + { + unsigned int eqIdx = (compIdx == wCompIdx) ? contiWEqIdx : contiNEqIdx; + storage[eqIdx] += volVars.molarDensity(phaseIdx) + * volVars.saturation(phaseIdx) + * volVars.moleFraction(phaseIdx, compIdx); + } + // this is only processed if one component mass balance equation + // is replaced by the total mass balance equation + if (replaceCompEqIdx < numComponents) + storage[replaceCompEqIdx] += + volVars.molarDensity(phaseIdx) + * volVars.saturation(phaseIdx); + } + storage *= volVars.porosity(); + } + else // mass-fraction formulation + { + for (unsigned int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) + { + for (unsigned int compIdx = contiCompIdx1_(); compIdx <= contiCompIdx2_(); ++compIdx) + { + unsigned int eqIdx = (compIdx == wCompIdx) ? contiWEqIdx : contiNEqIdx; + storage[eqIdx] += volVars.density(phaseIdx) + * volVars.saturation(phaseIdx) + * volVars.massFraction(phaseIdx, compIdx); + } + // this is only processed if one component mass balance equation + // is replaced by the total mass balance equation + if (replaceCompEqIdx < numComponents) + storage[replaceCompEqIdx] += + volVars.density(phaseIdx) + * volVars.saturation(phaseIdx); + } + storage *= volVars.porosity(); + } + } + + /*! + * \brief Evaluates the total flux of all conservation quantities + * over a face of a sub-control volume. + * + * \param flux The flux over the sub-control-volume face for each component + * \param fIdx The index of the sub-control-volume face + * \param onBoundary Evaluate flux at inner sub-control-volume face or on a boundary face + */ + void computeFlux(PrimaryVariables &flux, const int fIdx, bool onBoundary=false) const + { + // update the flux variables + FluxVariables fluxVars; + fluxVars.update(this->problem_(), + this->element_(), + this->fvGeometry_(), + fIdx, + this->curVolVars_(), + onBoundary); + + flux = 0; + + // add advective fluxes + computeAdvectiveFlux(flux, fluxVars); + + // add diffusive fluxes + // TODO: dumux-course-task + // only add diffusive fluxes if diffusion is enabled + } + + /*! + * \brief Evaluates the advective mass flux of all components over + * a face of a sub-control volume. + * + * \param flux The flux over the sub-control-volume face for each component + * \param fluxVars The flux variables at the current sub-control-volume face + */ + void computeAdvectiveFlux(PrimaryVariables &flux, const FluxVariables &fluxVars) const + { + //////// + // advective fluxes of all components in all phases + //////// + + if(useMoles) // mole-fraction formulation + { + for (unsigned int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) + { + // data attached to upstream and the downstream vertices + // of the current phase + const VolumeVariables &up = + this->curVolVars_(fluxVars.upstreamIdx(phaseIdx)); + const VolumeVariables &dn = + this->curVolVars_(fluxVars.downstreamIdx(phaseIdx)); + + for (unsigned int compIdx = contiCompIdx1_(); compIdx <= contiCompIdx2_(); ++compIdx) + { + unsigned int eqIdx = (compIdx == wCompIdx) ? contiWEqIdx : contiNEqIdx; + // add advective flux of current component in current + // phase + if (massUpwindWeight_ > 0.0) + // upstream vertex + flux[eqIdx] += + fluxVars.volumeFlux(phaseIdx) + * massUpwindWeight_ + * up.molarDensity(phaseIdx) + * up.moleFraction(phaseIdx, compIdx); + if (massUpwindWeight_ < 1.0) + // downstream vertex + flux[eqIdx] += + fluxVars.volumeFlux(phaseIdx) + * (1 - massUpwindWeight_) + * dn.molarDensity(phaseIdx) + * dn.moleFraction(phaseIdx, compIdx); + + Valgrind::CheckDefined(fluxVars.volumeFlux(phaseIdx)); + Valgrind::CheckDefined(up.molarDensity(phaseIdx)); + Valgrind::CheckDefined(up.moleFraction(phaseIdx, compIdx)); + Valgrind::CheckDefined(dn.molarDensity(phaseIdx)); + Valgrind::CheckDefined(dn.moleFraction(phaseIdx, compIdx)); + } + // flux of the total mass balance; + // this is only processed if one component mass balance equation + // is replaced by a total mass balance equation + if (replaceCompEqIdx < numComponents) + { + // upstream vertex + if (massUpwindWeight_ > 0.0) + flux[replaceCompEqIdx] += + fluxVars.volumeFlux(phaseIdx) + * massUpwindWeight_ + * up.molarDensity(phaseIdx); + // downstream vertex + if (massUpwindWeight_ < 1.0) + flux[replaceCompEqIdx] += + fluxVars.volumeFlux(phaseIdx) + * (1 - massUpwindWeight_) + * dn.molarDensity(phaseIdx); + Valgrind::CheckDefined(fluxVars.volumeFlux(phaseIdx)); + Valgrind::CheckDefined(up.molarDensity(phaseIdx)); + Valgrind::CheckDefined(dn.molarDensity(phaseIdx)); + + } + + } + } + else // mass-fraction formulation + { + for (unsigned int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) + { + // data attached to upstream and downstream vertices + // of the current phase + const VolumeVariables &up = + this->curVolVars_(fluxVars.upstreamIdx(phaseIdx)); + const VolumeVariables &dn = + this->curVolVars_(fluxVars.downstreamIdx(phaseIdx)); + + for (unsigned int compIdx = contiCompIdx1_(); compIdx <= contiCompIdx2_(); ++compIdx) + { + unsigned int eqIdx = (compIdx == wCompIdx) ? contiWEqIdx : contiNEqIdx; + // add advective flux of current component in current + // phase + if (massUpwindWeight_ > 0.0) + // upstream vertex + flux[eqIdx] += + fluxVars.volumeFlux(phaseIdx) + * massUpwindWeight_ + * up.density(phaseIdx) + * up.massFraction(phaseIdx, compIdx); + if (massUpwindWeight_ < 1.0) + // downstream vertex + flux[eqIdx] += + fluxVars.volumeFlux(phaseIdx) + * (1 - massUpwindWeight_) + * dn.density(phaseIdx) + * dn.massFraction(phaseIdx, compIdx); + + Valgrind::CheckDefined(fluxVars.volumeFlux(phaseIdx)); + Valgrind::CheckDefined(up.density(phaseIdx)); + Valgrind::CheckDefined(up.massFraction(phaseIdx, compIdx)); + Valgrind::CheckDefined(dn.density(phaseIdx)); + Valgrind::CheckDefined(dn.massFraction(phaseIdx, compIdx)); + } + // flux of the total mass balance; + // this is only processed if one component mass balance equation + // is replaced by a total mass balance equation + if (replaceCompEqIdx < numComponents) + { + // upstream vertex + if (massUpwindWeight_ > 0.0) + flux[replaceCompEqIdx] += + fluxVars.volumeFlux(phaseIdx) + * massUpwindWeight_ + * up.density(phaseIdx); + // downstream vertex + if (massUpwindWeight_ < 1.0) + flux[replaceCompEqIdx] += + fluxVars.volumeFlux(phaseIdx) + * (1 - massUpwindWeight_) + * dn.density(phaseIdx); + Valgrind::CheckDefined(fluxVars.volumeFlux(phaseIdx)); + Valgrind::CheckDefined(up.density(phaseIdx)); + Valgrind::CheckDefined(dn.density(phaseIdx)); + + } + + } + } + } + + /*! + * \brief Evaluates the diffusive mass flux of all components over + * a face of a sub-control volume. + * + * \param flux The flux over the sub-control-volume face for each component + * \param fluxVars The flux variables at the current sub-control-volume face + */ + void computeDiffusiveFlux(PrimaryVariables &flux, const FluxVariables &fluxVars) const + + { + if(useMoles) // mole-fraction formulation + { + // add diffusive flux of gas component in liquid phase + Scalar tmp = -(fluxVars.moleFractionGrad(wPhaseIdx)*fluxVars.face().normal) + * fluxVars.porousDiffCoeff(wPhaseIdx) + * fluxVars.molarDensity(wPhaseIdx); + // add the diffusive fluxes only to the component mass balance + if (replaceCompEqIdx != contiNEqIdx) + flux[contiNEqIdx] += tmp; + if (replaceCompEqIdx != contiWEqIdx) + flux[contiWEqIdx] -= tmp; + + // add diffusive flux of liquid component in non-wetting phase + tmp = -(fluxVars.moleFractionGrad(nPhaseIdx)*fluxVars.face().normal) + * fluxVars.porousDiffCoeff(nPhaseIdx) + * fluxVars.molarDensity(nPhaseIdx); + // add the diffusive fluxes only to the component mass balance + if (replaceCompEqIdx != contiWEqIdx) + flux[contiWEqIdx] += tmp; + if (replaceCompEqIdx != contiNEqIdx) + flux[contiNEqIdx] -= tmp; + } + else // mass-fraction formulation + { + // add diffusive flux of gas component in wetting phase + Scalar tmp = -(fluxVars.moleFractionGrad(wPhaseIdx)*fluxVars.face().normal) + * fluxVars.porousDiffCoeff(wPhaseIdx) + * fluxVars.molarDensity(wPhaseIdx); + // add the diffusive fluxes to the component mass balance and total mass balance + if (replaceCompEqIdx < numComponents) + { + flux[replaceCompEqIdx] += tmp * FluidSystem::molarMass(nCompIdx); + flux[replaceCompEqIdx] -= tmp * FluidSystem::molarMass(wCompIdx); + } + if (replaceCompEqIdx != contiWEqIdx) + { + flux[contiWEqIdx] -= tmp * FluidSystem::molarMass(wCompIdx); + } + if (replaceCompEqIdx != contiNEqIdx) + { + flux[contiNEqIdx] += tmp * FluidSystem::molarMass(nCompIdx); + } + + // add diffusive fluxes of liquid component in non-wetting phase + tmp = -(fluxVars.moleFractionGrad(nPhaseIdx)*fluxVars.face().normal) + * fluxVars.porousDiffCoeff(nPhaseIdx) + * fluxVars.molarDensity(nPhaseIdx); + // add the diffusive fluxes to the component mass balance and total mass balance + if (replaceCompEqIdx < numComponents) + { + flux[replaceCompEqIdx] += tmp * FluidSystem::molarMass(wCompIdx); + flux[replaceCompEqIdx] -= tmp * FluidSystem::molarMass(nCompIdx); + } + if (replaceCompEqIdx != contiWEqIdx) + { + flux[contiWEqIdx] += tmp * FluidSystem::molarMass(wCompIdx); + } + if (replaceCompEqIdx != contiNEqIdx) + { + flux[contiNEqIdx] -= tmp * FluidSystem::molarMass(nCompIdx); + } + } + } + + protected: + void evalPhaseStorage_(const int phaseIdx) + { + if(useMoles) // mole-fraction formulation + { + // evaluate the storage terms of a single phase + for (int i=0; i < this->fvGeometry_().numScv; i++) { + PrimaryVariables &storage = this->storageTerm_[i]; + const ElementVolumeVariables &elemVolVars = this->curVolVars_(); + const VolumeVariables &volVars = elemVolVars[i]; + + // compute storage term of all components within all phases + storage = 0; + for (int compIdx = 0; compIdx < numComponents; ++compIdx) + { + int eqIdx = (compIdx == wCompIdx) ? contiWEqIdx : contiNEqIdx; + storage[eqIdx] += volVars.molarDensity(phaseIdx) + * volVars.saturation(phaseIdx) + * volVars.moleFraction(phaseIdx, compIdx); + } + + storage *= volVars.porosity(); + storage *= this->fvGeometry_().subContVol[i].volume; + } + } + else // mass-fraction formulation + { + // evaluate the storage terms of a single phase + for (int i=0; i < this->fvGeometry_().numScv; i++) { + PrimaryVariables &storage = this->storageTerm_[i]; + const ElementVolumeVariables &elemVolVars = this->curVolVars_(); + const VolumeVariables &volVars = elemVolVars[i]; + + // compute storage term of all components within all phases + storage = 0; + for (int compIdx = 0; compIdx < numComponents; ++compIdx) + { + int eqIdx = (compIdx == wCompIdx) ? contiWEqIdx : contiNEqIdx; + storage[eqIdx] += volVars.density(phaseIdx) + * volVars.saturation(phaseIdx) + * volVars.massFraction(phaseIdx, compIdx); + } + + storage *= volVars.porosity(); + storage *= this->fvGeometry_().subContVol[i].volume; + } + } + } + + /*! + * \brief Returns the equation index of the first mass-balance equation + * of the component (used for loops) + * + * Returns the equation index of the first mass-balance equation + * of the component (used for loops) if one component mass balance + * is replaced by the total mass balance, this is the index + * of the remaining component mass-balance equation. + */ + unsigned int contiCompIdx1_() const { + switch (replaceCompEqIdx) + { + case contiWEqIdx: return contiNEqIdx; + case contiNEqIdx: return contiWEqIdx; + default: return 0; + } + } + + /*! + * \brief Returns the equation index of the second mass balance + * of the component (used for loops) + * + * Returns the equation index of the second mass balance + * of the component (used for loops) + * if one component mass balance is replaced by the total mass balance + * (replaceCompEqIdx < 2), this index is the same as contiCompIdx1(). + */ + unsigned int contiCompIdx2_() const { + switch (replaceCompEqIdx) + { + case contiWEqIdx: return contiNEqIdx; + case contiNEqIdx: return contiWEqIdx; + default: return numComponents-1; + } + } + + Implementation *asImp_() + { return static_cast (this); } + const Implementation *asImp_() const + { return static_cast (this); } + + private: + Scalar massUpwindWeight_; + + // TODO: dumux-course-task + // add private member enableDiffusion_ +}; + +} // end namespace Dumux + +#endif diff --git a/tutorial/ex2/mymateriallaw.hh b/tutorial/ex2/mymateriallaw.hh new file mode 100644 index 0000000000000000000000000000000000000000..662ee49fa3aac3ce3c240c18d4f6df1f473c547b --- /dev/null +++ b/tutorial/ex2/mymateriallaw.hh @@ -0,0 +1,115 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/***************************************************************************** + * See the file COPYING for full copying permissions. * + * * + * This program is free software: you can redistribute it and/or modify * + * it under the terms of the GNU General Public License as published by * + * the Free Software Foundation, either version 2 of the License, or * + * (at your option) any later version. * + * * + * This program is distributed in the hope that it will be useful, * + * but WITHOUT ANY WARRANTY; without even the implied warranty of * + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * + * GNU General Public License for more details. * + * * + * You should have received a copy of the GNU General Public License * + * along with this program. If not, see . * + *****************************************************************************/ +/*! + * \file + * + * \brief Implementation of the capillary pressure and + * relative permeability <-> saturation relations. + * + */ +#ifndef DUMUX_MY_MATERIAL_LAW_HH +#define DUMUX_MY_MATERIAL_LAW_HH + +#include +#include + +namespace Dumux +{ +/*! + * \ingroup fluidmatrixinteractionslaws + * \note a simple material law using the BrooksCoreyParams + */ +template > +class MyMaterialLaw +{ +public: + typedef ParamsT Params; + typedef typename Params::Scalar Scalar; + + /*! + * \brief The capillary pressure-saturation curve + * \param swe saturation of the wetting phase \f$\mathrm{[\overline{S}_w]}\f$ + * \param params A container object that is populated with the appropriate coefficients for the respective law. + * Therefore, in the (problem specific) spatialParameters first, the material law is chosen, + and then the params container is constructed accordingly. Afterwards the values are set there, too. + * \return capillary pressure + * TODO: dumux-course-task + * Implement the pc(swe) function + */ + static Scalar pc(const Params ¶ms, Scalar swe) + { + return 0.0; + } + + /*! + * \brief The relative permeability for the wetting phase of + * the medium implied by the Brooks-Corey + * parameterization. + * + * \param swe The mobile saturation of the wetting phase. + * \param params A container object that is populated with the appropriate coefficients for the respective law. + * Therefore, in the (problem specific) spatialParameters first, the material law is chosen, + * and then the params container is constructed accordingly. Afterwards the values are set there, too. + * \return Relative permeability of the wetting phase calculated as implied by Brooks & Corey. + * + * \note Instead of undefined behaviour if pc is not in the valid range, we return a valid number, + * by clamping the input. + */ + static Scalar krw(const Params ¶ms, Scalar swe) + { + using std::pow; + using std::min; + using std::max; + + swe = min(max(swe, 0.0), 1.0); // the equation below is only defined for 0.0 <= swe <= 1.0 + + return pow(swe, 2.0/params.lambda() + 3); + } + + /*! + * \brief The relative permeability for the non-wetting phase of + * the medium as implied by the Brooks-Corey + * parameterization. + * + * \param swe The mobile saturation of the wetting phase. + * \param params A container object that is populated with the appropriate coefficients for the respective law. + * Therefore, in the (problem specific) spatialParameters first, the material law is chosen, and then the params container + * is constructed accordingly. Afterwards the values are set there, too. + * \return Relative permeability of the non-wetting phase calculated as implied by Brooks & Corey. + * + * \note Instead of undefined behaviour if pc is not in the valid range, we return a valid number, + * by clamping the input. + */ + static Scalar krn(const Params ¶ms, Scalar swe) + { + using std::pow; + using std::min; + using std::max; + + swe = min(max(swe, 0.0), 1.0); // the equation below is only defined for 0.0 <= swe <= 1.0 + + const Scalar exponent = 2.0/params.lambda() + 1; + const Scalar tmp = 1.0 - swe; + return tmp*tmp*(1.0 - pow(swe, exponent)); + } +}; + +} // end namespace Dumux + +#endif diff --git a/tutorial/ex3/CMakeLists.txt b/tutorial/ex3/CMakeLists.txt new file mode 100644 index 0000000000000000000000000000000000000000..447a9a5f87611014c35c6e2d39a9a2cd45e45bb0 --- /dev/null +++ b/tutorial/ex3/CMakeLists.txt @@ -0,0 +1,17 @@ +# executables for exercise part a & b +dune_add_test(NAME ex3_a + SOURCES ex3_a.cc) +dune_add_test(NAME ex3_b + SOURCES ex3_b.cc) + +# add a symlink for the input file +dune_symlink_to_source_files(FILES "ex3_a.input" "ex3_b.input") + +#install sources +install(FILES +ex3_a.cc +ex3_a_problem.hh +ex3_b.cc +ex3_b_problem.hh +spatialparams.hh +DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dumux/tutorial/ex3) diff --git a/tutorial/ex3/README.md b/tutorial/ex3/README.md new file mode 100644 index 0000000000000000000000000000000000000000..e803a7e83855abeb0070d9555e19e118ed8f5f5b --- /dev/null +++ b/tutorial/ex3/README.md @@ -0,0 +1,227 @@ +# Exercise #3 (DuMuX course) + +The aim of this exercise is to get familiar with the _DuMuX_ way of implementing new components (fluids) and fluid systems (mixtures). In the scope of this exercise, a new fictitious component is implemented (exercise _3a_) as well as its mixture with water (exercise _3b_). + +## Problem set-up + +The domain has a size of 60 x 60 m and contains two low-permeable lenses. Initially, the domain is fully water saturated and the fictitious component is injected through the middle portion of the upper boundary by means of a Neumann boundary condition. The remaining parts of the upper and the entire lower boundary are Neumann no-flow while on the two lateral sides Dirichlet boundary conditions are applied (hydrostatic conditions for the pressure and zero saturation). + +![](../extradoc/exercise3_setup.png) + + +## Preparing the exercise + +* Navigate to the directory `dumux/tutorial/ex3` + +### 1. Getting familiar with the code + +Locate all the files you will need for this exercise +* The __main file__ for part a: `ex3_a.cc` +* The __input file__ for part a: `ex3_a.input` +* The __problem file__ for part a: `ex3_a_problem.hh` +* The __main file__ for part b: `ex3_b.cc` +* The __input file__ for part b: `ex3_b.input` +* The __problem file__ for part b: `ex3_b_problem.hh` +* The __spatial parameters file__: `spatialparams.hh` + +Furthermore you will find the following folders: +* `binarycoefficients`: Stores headers containing data/methods on binary mixtures +* `components`: Stores headers containing data/methods on pure components +* `fluidsystems`: Stores headers containing data/methods on mixtures of pure components. Uses methods from `binarycoefficients`. + +To see more components, fluidsystems and binarycoefficients implementations, have a look at the folder `dumux/material`. + +### 2. Implement a new component + +In the following the basic steps required to set the desired fluid system are outlined. Here, this is done in the __problem file__, i.e. for this part of the exercise the code shown below is taken from the `ex3_a_problem.hh` file. + +In this part of the exercise we will consider a system consisting of two immiscible phases. Therefore, the _TypeTag_ for this problem derives from the _BoxTwoP_ _TypeTag_ (for a cell-centered scheme, you would derive from _CCTwoP_). In order to be able to derive from this _TypeTag_, The declaration of the _BoxTwoP_ _TypeTag_ is found in the file `dumux/porousmediumflow/2p/implicit/model.hh`, which has been included in line 28: + +```c++ +// The numerical model +#include +``` + + Additionally we derive from the _ExerciseThreeSpatialParams_ _TypeTag_ in order to set all the types related to the spatial parameters also for our new _TypeTag_ _ExerciseThreeProblem_. In this case, only one property is set for _ExerciseThreeSpatialParams_ (see lines 43 to 60 in `spatialparams.hh`). Alternatively, we could have defined this also here in the problem without defining a new type tag for the spatial parameters. However, in case you want to define several properties related to the spatial parameters it is good practice to define a separate _TypeTag_ and derive from this, as it is done here. + +```c++ +NEW_TYPE_TAG(ExerciseThreeProblem, INHERITS_FROM(BoxTwoP, ExerciseThreeSpatialParams)); +``` + +As wetting phase we want to use water and we want to precompute tables on which the properties are then interpolated in order to save computational time. Thus, in a first step we have to include the following headers: + +```c++ +// The water component +#include +#include +``` +The non-wetting phase will be our new component, where we want to implement an incompressible and a compressible variant. The respective headers are prepared, but still incomplete. The compressible variant is still commented so that compilation does not fail when finishing the incompressible variant. + +```c++ +// The components that will be created in this exercise +#include "components/myincompressiblecomponent.hh" +// #include "components/mycompressiblecomponent.hh" +``` +As mentioned above, we want to simulate two non-mixing components. The respective fluid system is found in: + +```c++ +// The two-phase immiscible fluid system +#include +``` + +This fluid system expects __phases__ as input and so far we have only included the components, which contain data on the pure component for all physical states. Thus, we need to include + +```c++ +// We will only have liquid phases Here +#include +``` + +which creates a _liquid phase_ from a given component. Finally, using all of the included classes we set the the fluid system property by choosing a liquid phase consisting of the incompressible fictitious component as non-wetting phase and tabulated water as the wetting phase in the immiscible fluid system: + + +```c++ +// we use the immiscible fluid system here +SET_PROP(ExerciseThreeProblem, FluidSystem) +{ +private: + typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar; + typedef TabulatedComponent> TabulatedH2O; + typedef typename FluidSystems::LiquidPhase WettingPhase; + /*! + * Uncomment first line and comment second line for using the incompressible component + * Uncomment second line and comment first line for using the compressible component + */ + typedef typename FluidSystems::LiquidPhase > NonWettingPhase; + // typedef typename FluidSystems::LiquidPhase > NonWettingPhase; + +public: + typedef typename FluidSystems::TwoPImmiscible type; +}; +``` + +### 2.1. Incompressible component + +Open the file `myincompressiblecomponent.hh`. You can see in line 40 that a component should always derive from the _Component_ class (see `dumux/material/components/component.hh`), which defines the interface of a _DuMuX_ component with all possibly required functions to be overloaded by the actual implementation. + +```c++ +/*! + * \ingroup Components + * \brief A ficitious component to be implemented in exercise 3. + * + * \tparam Scalar The type used for scalar values + */ +template +class MyIncompressibleComponent : public Component > +``` + +__Task__: + +Implement an incompressible component into the file `myincompressiblecomponent.hh`, which has the following specifications: + +| Parameter | unit | value | +| -----| --------| -------- | +| $`M`$ | $`Kg/mol`$ | $`131.39 \cdot 10^{-3}`$ | +| $`\rho_{liquid}`$ | $`Kg/m^3`$ | $`1460`$ | +| $`\mu_{liquid}`$ | $`Pa \cdot s`$ | $`5.7 \cdot 10^{-4}`$ | + +In order to do so, have a look at the file `dumux/material/components/component.hh` to see how the interfaces are defined and overload them accordingly. + +Execute the program by changing to the build-directory + +```bash +cd build-cmake/tutorial/ex3 +``` + +and by compiling the executable for part a + +```bash +make ex3_a +``` + +and by typing + +```bash +./ex3_a ex3_a.input +``` + +The saturation distribution at the final simulation time should look like this: + +![](../extradoc/exercise3_a_solution.png) + +### 2.2. Compressible component + +We now want to implement a pressure-dependent density for our component. Open the file `mycompressiblecomponent.hh` and copy in the functions you implemented for the incompressible variant. Now substitute the method that returns the density by the following expression: + +$`\rho_{MyComp} = \rho_{min} + \frac{ \rho_{max} - \rho_{min} }{ 1 + \rho_{min}*e^{-1.0*k*(\rho_{max} - \rho_{min})*p} } `$, + +where $`p`$ is the pressure and $`\rho_{min} = 1440 `$, $`\rho_{max} = 1480 `$ and $`k = 5 \cdot 10^{-7} `$. Also, make sure the header is included in the `ex3_a_problem.hh` file by uncommenting line 42. Furthermore, the new component has to be set as the non-wetting phase in the fluid system, i.e. comment line 81 and uncomment line 82. The non-wetting density distribution at the final simulation time should look like this: + +![](../extradoc/exercise3_a_solution2.png) + +### 3. Implement a new fluid system + +The problem file for this part of the exercise is `ex3_b_problem.hh`. We now want to implement a new fluid system consisting of two liquid phases, which are water and the previously implemented compressible component. We will consider compositional effects, which is why we now have to derive our _TypeTag_ from the _BoxTwoPTwoC_ _TypeTag_: + +```c++ +// The numerical model +#include +``` + +```c++ +// Create a new type tag for the problem +NEW_TYPE_TAG(ExerciseThreeProblem, INHERITS_FROM(BoxTwoPTwoC, ExerciseThreeSpatialParams)); +``` + +The new fluid system is to be implemented in the file `fluidsystems/h2omycompressiblecomponent.hh`. This is already included in the problem and the fluid system property is set accordingly. + +```c++ +// The fluid system that is created in this exercise +#include "fluidsystems/h2omycompressiblecomponent.hh" +``` + +```c++ +// The fluid system property +SET_PROP(ExerciseThreeProblem, FluidSystem) +{ +private: + typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar; +public: + typedef FluidSystems::H2OMyCompressibleComponent type; +}; +``` + +In the `fluidsystems/h2omycompressiblecomponent.hh` file, your implemented component and the binary coefficient files are already included. + +```c++ +// the ficitious component that was created in exercise 3a +#include + +// the binary coefficients corresponding to this fluid system +#include +``` + +__Task__: + +Under the assumption that one molecule of _MyCompressibleComponent_ displaces exactly one molecule of water, the water phase density can be expressed as follows: + +$` \rho_{w} = \frac{ \rho_{w, pure} }{ M_{H_2O} }*(M_{H_2O}*x_{H_2O} + M_{MyComponent}*x_{MyComponent}) `$ + +Implement this dependency in the _density()_ method in the fluid system. Then, execute the program by changing to the build-directory + +```bash +cd build-cmake/tutorial/ex3 +``` + +and by compiling the executable for part b + +```bash +make ex3_b +``` + +and by typing + +```bash +./ex3_b ex3_b.input +``` + +You will observe an error message and an abortion of the program. This is due to the fact that in order for the constraint solver and other mechanisms in the two-phase two-component model to work, two additional functionalities in the component have to be implemented. The model has to know whether or not the liquid pure component is compressible and it needs the vapour pressure. As in the previous exercise, check the `dumux/material/components/component.hh` file for these two functions and implement them into `mycompressiblecomponent.hh`. For the vapour pressure, use a value of $`3900`$ Pa. diff --git a/tutorial/ex3/binarycoefficients/h2omycompressiblecomponent.hh b/tutorial/ex3/binarycoefficients/h2omycompressiblecomponent.hh new file mode 100644 index 0000000000000000000000000000000000000000..843a69a675bc1b8095a850c94f07a8ca65232228 --- /dev/null +++ b/tutorial/ex3/binarycoefficients/h2omycompressiblecomponent.hh @@ -0,0 +1,74 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/***************************************************************************** + * See the file COPYING for full copying permissions. * + * * + * This program is free software: you can redistribute it and/or modify * + * it under the terms of the GNU General Public License as published by * + * the Free Software Foundation, either version 2 of the License, or * + * (at your option) any later version. * + * * + * This program is distributed in the hope that it will be useful, * + * but WITHOUT ANY WARRANTY; without even the implied warranty of * + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * + * GNU General Public License for more details. * + * * + * You should have received a copy of the GNU General Public License * + * along with this program. If not, see . * + *****************************************************************************/ +/*! + * \file + * + * \brief Binary coefficients for water and a ficticious component implemented in tutorial exercise 3a. + */ +#ifndef DUMUX_BINARY_COEFF_H2O_MYCOMPRESSIBLECOMPONENT_HH +#define DUMUX_BINARY_COEFF_H2O_MYCOMPRESSIBLECOMPONENT_HH + +namespace Dumux +{ +namespace BinaryCoeff +{ + +/*! + * \brief Binary coefficients for water and a ficticious component implemented in tutorial exercise 3a + * The implementation of the missing methods in this file is part of exercise 3b. + */ +class H2O_MyCompressibleComponent +{ +public: + /*! + * \brief Henry coefficient \f$[N/m^2]\f$ for the fictitous component in liquid water. + */ + template + static Scalar henryMyCompressibleComponentInWater(Scalar temperature) + { + Scalar dumuxH = 1.5e-1 / 101.325; // unit [(mol/m^3)/Pa] + dumuxH *= 18.02e-6; //multiplied by molar volume of reference phase = water + return 1.0/dumuxH; // [Pa] + } + + /*! + * \brief Henry coefficient \f$[N/m^2]\f$ for water in the ficticious component. + */ + template + static Scalar henryWaterInMyCompressibleComponent(Scalar temperature) + { + // arbitrary + return 1.0e8; // [Pa] + } + + /*! + * \brief Diffusion coefficient [m^2/s] for my ficticious component in liquid water or vice versa. + */ + template + static Scalar liquidDiffCoeff(Scalar temperature, Scalar pressure) + { + // arbitrary + return 1.e-9; + } +}; + +} +} // end namespace + +#endif diff --git a/tutorial/ex3/components/mycompressiblecomponent.hh b/tutorial/ex3/components/mycompressiblecomponent.hh new file mode 100644 index 0000000000000000000000000000000000000000..65d48fe25a68c0cabba02fa79419627f540c262c --- /dev/null +++ b/tutorial/ex3/components/mycompressiblecomponent.hh @@ -0,0 +1,58 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/***************************************************************************** + * See the file COPYING for full copying permissions. * + * * + * This program is free software: you can redistribute it and/or modify * + * it under the terms of the GNU General Public License as published by * + * the Free Software Foundation, either version 2 of the License, or * + * (at your option) any later version. * + * * + * This program is distributed in the hope that it will be useful, * + * but WITHOUT ANY WARRANTY; without even the implied warranty of * + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * + * GNU General Public License for more details. * + * * + * You should have received a copy of the GNU General Public License * + * along with this program. If not, see . * + *****************************************************************************/ +/*! + * \file + * \ingroup Components + * \brief A ficitious component to be implemented in exercise 3. + */ +#ifndef DUMUX_MYCOMPRESSIBLECOMPONENT_HH +#define DUMUX_MYCOMPRESSIBLECOMPONENT_HH + +#include +#include + + +namespace Dumux +{ +/*! + * \ingroup Components + * \brief A ficitious component to be implemented in exercise 3. + * + * \tparam Scalar The type used for scalar values + */ +template +class MyCompressibleComponent : public Component > +{ + +public: + /*! + * \brief A human readable name for MyCompressibleComponent. + */ + static std::string name() + { return "MyCompressibleComponent"; } + + /*! + * TODO: Copy the methods implemented in MyIncompressibleComponent and substitute + * the density calculation by the expression given in the exercise description. + */ +}; + +} // end namespace + +#endif diff --git a/tutorial/ex3/components/myincompressiblecomponent.hh b/tutorial/ex3/components/myincompressiblecomponent.hh new file mode 100644 index 0000000000000000000000000000000000000000..9809100b21d95e0ae207858f4d61cfe95c4cdce3 --- /dev/null +++ b/tutorial/ex3/components/myincompressiblecomponent.hh @@ -0,0 +1,56 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/***************************************************************************** + * See the file COPYING for full copying permissions. * + * * + * This program is free software: you can redistribute it and/or modify * + * it under the terms of the GNU General Public License as published by * + * the Free Software Foundation, either version 2 of the License, or * + * (at your option) any later version. * + * * + * This program is distributed in the hope that it will be useful, * + * but WITHOUT ANY WARRANTY; without even the implied warranty of * + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * + * GNU General Public License for more details. * + * * + * You should have received a copy of the GNU General Public License * + * along with this program. If not, see . * + *****************************************************************************/ +/*! + * \file + * \ingroup Components + * \brief A ficitious component to be implemented in tutorial exercise 3. + */ +#ifndef DUMUX_MYINCOMPRESSIBLECOMPONENT_HH +#define DUMUX_MYINCOMPRESSIBLECOMPONENT_HH + +#include +#include + + +namespace Dumux +{ +/*! + * \ingroup Components + * \brief A ficitious component to be implemented in exercise 3. + * + * \tparam Scalar The type used for scalar values + */ +template +class MyIncompressibleComponent : public Component > +{ +public: + /*! + * \brief A human readable name for MyIncompressibleComponent. + */ + static std::string name() + { return "MyIncompressibleComponent"; } + + /*! + * TODO: Implement the methods for the component data given in the exercise description. + */ +}; + +} // end namespace + +#endif diff --git a/tutorial/ex3/components/plotdensityfunction.py b/tutorial/ex3/components/plotdensityfunction.py new file mode 100644 index 0000000000000000000000000000000000000000..2b256689339b8c094d4bd3afb235be2b57fbb214 --- /dev/null +++ b/tutorial/ex3/components/plotdensityfunction.py @@ -0,0 +1,19 @@ +#!usr/bin/env python +import numpy as np +import matplotlib.pyplot as plt + +# function to calculate rho dependent on pressure +rho_min = 1440; +rho_max = 1480; +k = 5e-7; + +def rho(p): + return rho_min + (rho_max - rho_min)/(1 + rho_min*np.exp(-1.0*k*(rho_max - rho_min)*p)); + +# sample pressure in range (1e4, 1e7) and compute corresponding densities +p = np.logspace(4, 7, 100) +r = rho(p) + +# plot density vs. pressure +plt.semilogx(p, r) +plt.show() diff --git a/tutorial/ex3/ex3_a.cc b/tutorial/ex3/ex3_a.cc new file mode 100644 index 0000000000000000000000000000000000000000..1761630a24e639fcace076f2b205bf70fb713cce --- /dev/null +++ b/tutorial/ex3/ex3_a.cc @@ -0,0 +1,49 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/***************************************************************************** + * See the file COPYING for full copying permissions. * + * * + * This program is free software: you can redistribute it and/or modify * + * it under the terms of the GNU General Public License as published by * + * the Free Software Foundation, either version 2 of the License, or * + * (at your option) any later version. * + * * + * This program is distributed in the hope that it will be useful, * + * but WITHOUT ANY WARRANTY; without even the implied warranty of * + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * + * GNU General Public License for more details. * + * * + * You should have received a copy of the GNU General Public License * + * along with this program. If not, see . * + *****************************************************************************/ +/*! + * \file + * + * \brief Main file of exercise 3. + */ +#include +#include "ex3_a_problem.hh" +#include + +//! Prints a usage/help message if something goes wrong or the user asks for help +void usage(const char *progName, const std::string &errorMsg) +{ + std::cout + << "\nUsage: " << progName << " [options]\n"; + if (errorMsg.size() > 0) + std::cout << errorMsg << "\n"; + std::cout + << "\n" + << "The list of mandatory arguments for this program is:\n" + << "\t-TEnd The end of the simulation [s]\n" + << "\t-DtInitial The initial timestep size [s]\n" + << "\t-Grid.UpperRight The x-/y-coordinates of the grid's upper-right corner [m]\n" + << "\t-Grid.Cells The grid's x-/y-resolution\n" + << "\n"; +} + +int main(int argc, char** argv) +{ + typedef TTAG(ExerciseThreeProblem) TypeTag; + return Dumux::start(argc, argv, usage); +} diff --git a/tutorial/ex3/ex3_a.input b/tutorial/ex3/ex3_a.input new file mode 100644 index 0000000000000000000000000000000000000000..ce45509020c21dbf243fda9bb4b9522cb46b60d8 --- /dev/null +++ b/tutorial/ex3/ex3_a.input @@ -0,0 +1,10 @@ +[TimeManager] +TEnd = 20000 # duration of the simulation [s] +DtInitial = 10 # initial time step size [s] + +[Problem] +Name = exercise3_a + +[Grid] +UpperRight = 60 60 # x-/y-coordinates of the upper-right corner of the grid [m] +Cells = 60 60 # x-/y-resolution of the grid diff --git a/tutorial/ex3/ex3_a_problem.hh b/tutorial/ex3/ex3_a_problem.hh new file mode 100644 index 0000000000000000000000000000000000000000..897d5c57ca1b62d6b7d34009d274636aef88938e --- /dev/null +++ b/tutorial/ex3/ex3_a_problem.hh @@ -0,0 +1,276 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/***************************************************************************** + * See the file COPYING for full copying permissions. * + * * + * This program is free software: you can redistribute it and/or modify * + * it under the terms of the GNU General Public License as published by * + * the Free Software Foundation, either version 2 of the License, or * + * (at your option) any later version. * + * * + * This program is distributed in the hope that it will be useful, * + * but WITHOUT ANY WARRANTY; without even the implied warranty of * + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * + * GNU General Public License for more details. * + * * + * You should have received a copy of the GNU General Public License * + * along with this program. If not, see . * + *****************************************************************************/ +/*! + * \file + * + * \brief Tutorial problem for a fully coupled twophase box model. + */ +#ifndef DUMUX_EXERCISE_THREE_PROBLEM_HH +#define DUMUX_EXERCISE_THREE_PROBLEM_HH + +// The numerical model +#include + +// The base porous media box problem +#include + +// Spatially dependent parameters +#include "spatialparams.hh" + +// The water component +#include +#include + +// The components that will be created in this exercise +#include "components/myincompressiblecomponent.hh" +// #include "components/mycompressiblecomponent.hh" + +// We will only have liquid phases Here +#include + +// The two-phase immiscible fluid system +#include + +namespace Dumux{ +// Forward declaration of the problem class +template class ExerciseThreeProblem; + +namespace Properties { +// Create a new type tag for the problem +NEW_TYPE_TAG(ExerciseThreeProblem, INHERITS_FROM(BoxTwoP, ExerciseThreeSpatialParams)); + +// Set the "Problem" property +SET_TYPE_PROP(ExerciseThreeProblem, Problem, ExerciseThreeProblem); + +// Set grid and the grid creator to be used +#if HAVE_DUNE_ALUGRID +SET_TYPE_PROP(ExerciseThreeProblem, Grid, Dune::ALUGrid); +#elif HAVE_UG +SET_TYPE_PROP(ExerciseThreeProblem, Grid, Dune::UGGrid<2>); +#else +SET_TYPE_PROP(ExerciseThreeProblem, Grid, Dune::YaspGrid<2>); +#endif // HAVE_DUNE_ALUGRID + +// we use the immiscible fluid system here +SET_PROP(ExerciseThreeProblem, FluidSystem) +{ +private: + typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar; + typedef TabulatedComponent> TabulatedH2O; + typedef typename FluidSystems::LiquidPhase WettingPhase; + /*! + * Uncomment first line and comment second line for using the incompressible component + * Uncomment second line and comment first line for using the compressible component + */ + typedef typename FluidSystems::LiquidPhase > NonWettingPhase; + // typedef typename FluidSystems::LiquidPhase > NonWettingPhase; + +public: + typedef typename FluidSystems::TwoPImmiscible type; +}; + +// Disable gravity +SET_BOOL_PROP(ExerciseThreeProblem, ProblemEnableGravity, true); +} + +/*! + * \ingroup TwoPBoxModel + * + * \brief Tutorial problem for a fully coupled twophase box model. + */ +template +class ExerciseThreeProblem : public ImplicitPorousMediaProblem +{ + typedef ImplicitPorousMediaProblem ParentType; + typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar; + typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView; + + // Grid dimension + enum { dim = GridView::dimension, + dimWorld = GridView::dimensionworld + }; + + // Types from DUNE-Grid + typedef typename GridView::template Codim<0>::Entity Element; + typedef typename GridView::template Codim::Entity Vertex; + typedef typename GridView::Intersection Intersection; + typedef Dune::FieldVector GlobalPosition; + + // Dumux specific types + typedef typename GET_PROP_TYPE(TypeTag, TimeManager) TimeManager; + typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices; + typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables; + typedef typename GET_PROP_TYPE(TypeTag, BoundaryTypes) BoundaryTypes; + typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry; + +public: + ExerciseThreeProblem(TimeManager &timeManager, + const GridView &gridView) + : ParentType(timeManager, gridView) + , eps_(3e-6) + { +#if !(HAVE_DUNE_ALUGRID || HAVE_UG) + std::cout << "If you want to use simplices instead of cubes, install and use dune-ALUGrid or UGGrid." << std::endl; +#endif // !(HAVE_DUNE_ALUGRID || HAVE_UG) + + // initialize the tables for the water properties + std::cout << "Initializing the tables for the water properties" << std::endl; + TabulatedComponent>::init(/*tempMin=*/273.15, + /*tempMax=*/623.15, + /*numTemp=*/100, + /*pMin=*/0.0, + /*pMax=*/20e6, + /*numP=*/200); + + // set the depth of the bottom of the reservoir + depthBOR_ = this->bBoxMax()[dimWorld-1]; + } + + //! Specifies the problem name. This is used as a prefix for files + //! generated by the simulation. + std::string name() const + { + static const std::string name = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, std::string, Problem, Name); + return name; + } + + //! Returns true if a restart file should be written. + bool shouldWriteRestartFile() const + { return false; } + + //! Returns true if the current solution should be written to disk + //! as a VTK file + bool shouldWriteOutput() const + { + return this->timeManager().timeStepIndex() > 0 && + this->timeManager().timeStepIndex() % 1 == 0; + } + + //! Returns the temperature within a finite volume. We use constant + //! 10 degrees Celsius. + Scalar temperature() const + { return 283.15; } + + //! Specifies which kind of boundary condition should be used for + //! which equation for a finite volume on the boundary. + void boundaryTypes(BoundaryTypes &bcTypes, const Vertex &vertex) const + { + const GlobalPosition &globalPos = vertex.geometry().center(); + + // Dirichlet conditions on left & right boundary + if (globalPos[0] < eps_ || globalPos[0] > this->bBoxMax()[0] - eps_) + bcTypes.setAllDirichlet(); + else // neuman for the remaining boundaries + bcTypes.setAllNeumann(); + + } + + //! Evaluates the Dirichlet boundary conditions for a finite volume + //! on the grid boundary. Here, the 'values' parameter stores + //! primary variables. + void dirichlet(PrimaryVariables &values, const Vertex &vertex) const + { + const auto globalPos = vertex.geometry().center(); + values[Indices::pwIdx] = 200.0e3 + 9.81*1000*(depthBOR_ - globalPos[dimWorld-1]); // 200 kPa = 2 bar + values[Indices::snIdx] = 0.0; // 0 % oil saturation on left boundary + } + + //! Evaluates the boundary conditions for a Neumann boundary + //! segment. Here, the 'values' parameter stores the mass flux in + //! [kg/(m^2 * s)] in normal direction of each phase. Negative + //! values mean influx. + void neumann(PrimaryVariables &values, + const Element &element, + const FVElementGeometry &fvGeometry, + const Intersection &intersection, + int scvIdx, + int boundaryFaceIdx) const + { + const GlobalPosition &globalPos = + fvGeometry.boundaryFace[boundaryFaceIdx].ipGlobal; + Scalar up = this->bBoxMax()[dimWorld-1]; + // extraction of oil on the right boundary for approx. 1.e6 seconds + if (globalPos[dimWorld-1] > up - eps_ && globalPos[0] > 20 && globalPos[0] < 40) { + // oil outflux of 30 g/(m * s) on the right boundary. + values[Indices::contiWEqIdx] = 0; + values[Indices::contiNEqIdx] = -3e-2; + } else { + // no-flow on the remaining Neumann-boundaries. + values[Indices::contiWEqIdx] = 0; + values[Indices::contiNEqIdx] = 0; + } + } + + //! Evaluates the initial value for a control volume. For this + //! method, the 'values' parameter stores primary variables. + void initial(PrimaryVariables &values, + const Element &element, + const FVElementGeometry &fvGeometry, + int scvIdx) const + { + const auto globalPos = fvGeometry.subContVol[scvIdx].global; + values[Indices::pwIdx] = 200.0e3 + 9.81*1000*(depthBOR_ - globalPos[dimWorld-1]); // 200 kPa = 2 bar + values[Indices::snIdx] = 0.0; + } + + //! Evaluates the source term for all phases within a given + //! sub-control-volume. In this case, the 'values' parameter + //! stores the rate mass generated or annihilated per volume unit + //! in [kg / (m^3 * s)]. Positive values mean that mass is created. + void source(PrimaryVariables &values, + const Element &element, + const FVElementGeometry &fvGeometry, + int scvIdx) const + { + values[Indices::contiWEqIdx] = 0.0; + values[Indices::contiNEqIdx]= 0.0; + } + + /*! + * \brief Append all quantities of interest which can be derived + * from the solution of the current time step to the VTK + * writer. Adjust this in case of anisotropic permeabilities. + */ + void addOutputVtkFields() + { + // get the number of elements in the grid + unsigned numElements = this->gridView().size(/*codim=*/0); + + // create the scalar field required for the output of the lenses + typedef Dune::BlockVector > ScalarField; + ScalarField *isInLens = this->resultWriter().allocateManagedBuffer(numElements); + + // add data to the scalar field + for (const auto& element : elements(this->gridView())) + (*isInLens)[this->model().elementMapper().index(element)] = this->spatialParams().isInLens(element.geometry().center()); + + // attach to writer + this->resultWriter().attachCellData(*isInLens, "isInLens"); + } + +private: + // small epsilon value + Scalar eps_; + + // depth at the bottom of the reservoir + Scalar depthBOR_; +}; +} + +#endif diff --git a/tutorial/ex3/ex3_b.cc b/tutorial/ex3/ex3_b.cc new file mode 100644 index 0000000000000000000000000000000000000000..1a548af8123a57d03c26a745d23bb67889517bc7 --- /dev/null +++ b/tutorial/ex3/ex3_b.cc @@ -0,0 +1,49 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/***************************************************************************** + * See the file COPYING for full copying permissions. * + * * + * This program is free software: you can redistribute it and/or modify * + * it under the terms of the GNU General Public License as published by * + * the Free Software Foundation, either version 2 of the License, or * + * (at your option) any later version. * + * * + * This program is distributed in the hope that it will be useful, * + * but WITHOUT ANY WARRANTY; without even the implied warranty of * + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * + * GNU General Public License for more details. * + * * + * You should have received a copy of the GNU General Public License * + * along with this program. If not, see . * + *****************************************************************************/ +/*! + * \file + * + * \brief Main file of exercise 3. + */ +#include +#include "ex3_b_problem.hh" +#include + +//! Prints a usage/help message if something goes wrong or the user asks for help +void usage(const char *progName, const std::string &errorMsg) +{ + std::cout + << "\nUsage: " << progName << " [options]\n"; + if (errorMsg.size() > 0) + std::cout << errorMsg << "\n"; + std::cout + << "\n" + << "The list of mandatory arguments for this program is:\n" + << "\t-TEnd The end of the simulation [s]\n" + << "\t-DtInitial The initial timestep size [s]\n" + << "\t-Grid.UpperRight The x-/y-coordinates of the grid's upper-right corner [m]\n" + << "\t-Grid.Cells The grid's x-/y-resolution\n" + << "\n"; +} + +int main(int argc, char** argv) +{ + typedef TTAG(ExerciseThreeProblem) TypeTag; + return Dumux::start(argc, argv, usage); +} diff --git a/tutorial/ex3/ex3_b.input b/tutorial/ex3/ex3_b.input new file mode 100644 index 0000000000000000000000000000000000000000..16c1b6a5ece6e615f594c5ce4375ae3afd159b4d --- /dev/null +++ b/tutorial/ex3/ex3_b.input @@ -0,0 +1,10 @@ +[TimeManager] +TEnd = 20000 # duration of the simulation [s] +DtInitial = 10 # initial time step size [s] + +[Problem] +Name = exercise3_b + +[Grid] +UpperRight = 60 60 # x-/y-coordinates of the upper-right corner of the grid [m] +Cells = 60 60 # x-/y-resolution of the grid diff --git a/tutorial/ex3/ex3_b_problem.hh b/tutorial/ex3/ex3_b_problem.hh new file mode 100644 index 0000000000000000000000000000000000000000..766396f329f1344bfeffb8c81f24846cc6f6d1cf --- /dev/null +++ b/tutorial/ex3/ex3_b_problem.hh @@ -0,0 +1,256 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/***************************************************************************** + * See the file COPYING for full copying permissions. * + * * + * This program is free software: you can redistribute it and/or modify * + * it under the terms of the GNU General Public License as published by * + * the Free Software Foundation, either version 2 of the License, or * + * (at your option) any later version. * + * * + * This program is distributed in the hope that it will be useful, * + * but WITHOUT ANY WARRANTY; without even the implied warranty of * + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * + * GNU General Public License for more details. * + * * + * You should have received a copy of the GNU General Public License * + * along with this program. If not, see . * + *****************************************************************************/ +/*! + * \file + * + * \brief Tutorial problem for a fully coupled twophase box model. + */ +#ifndef DUMUX_EXERCISE_THREE_PROBLEM_HH +#define DUMUX_EXERCISE_THREE_PROBLEM_HH + +// The numerical model +#include + +// The base porous media box problem +#include + +// Spatially dependent parameters +#include "spatialparams.hh" + +// The fluid system that is created in this exercise +#include "fluidsystems/h2omycompressiblecomponent.hh" + +namespace Dumux{ +// Forward declaration of the problem class +template class ExerciseThreeProblem; + +namespace Properties { +// Create a new type tag for the problem +NEW_TYPE_TAG(ExerciseThreeProblem, INHERITS_FROM(BoxTwoPTwoC, ExerciseThreeSpatialParams)); + +// Set the "Problem" property +SET_TYPE_PROP(ExerciseThreeProblem, Problem, ExerciseThreeProblem); + +// Set grid and the grid creator to be used +#if HAVE_DUNE_ALUGRID +SET_TYPE_PROP(ExerciseThreeProblem, Grid, Dune::ALUGrid); +#elif HAVE_UG +SET_TYPE_PROP(ExerciseThreeProblem, Grid, Dune::UGGrid<2>); +#else +SET_TYPE_PROP(ExerciseThreeProblem, Grid, Dune::YaspGrid<2>); +#endif // HAVE_DUNE_ALUGRID + + // The fluid system property +SET_PROP(ExerciseThreeProblem, FluidSystem) +{ +private: + typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar; +public: + typedef FluidSystems::H2OMyCompressibleComponent type; +}; + +// Disable gravity +SET_BOOL_PROP(ExerciseThreeProblem, ProblemEnableGravity, true); +} + +/*! + * \ingroup TwoPBoxModel + * + * \brief Tutorial problem for a fully coupled twophase box model. + */ +template +class ExerciseThreeProblem : public ImplicitPorousMediaProblem +{ + typedef ImplicitPorousMediaProblem ParentType; + typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar; + typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView; + + // Grid dimension + enum { dim = GridView::dimension, + dimWorld = GridView::dimensionworld + }; + + // Types from DUNE-Grid + typedef typename GridView::template Codim<0>::Entity Element; + typedef typename GridView::template Codim::Entity Vertex; + typedef typename GridView::Intersection Intersection; + typedef Dune::FieldVector GlobalPosition; + + // Dumux specific types + typedef typename GET_PROP_TYPE(TypeTag, TimeManager) TimeManager; + typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices; + typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables; + typedef typename GET_PROP_TYPE(TypeTag, BoundaryTypes) BoundaryTypes; + typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry; + typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem; + +public: + ExerciseThreeProblem(TimeManager &timeManager, + const GridView &gridView) + : ParentType(timeManager, gridView) + , eps_(3e-6) + { +#if !(HAVE_DUNE_ALUGRID || HAVE_UG) + std::cout << "If you want to use simplices instead of cubes, install and use dune-ALUGrid or UGGrid." << std::endl; +#endif // !(HAVE_DUNE_ALUGRID || HAVE_UG) + + // initialize the fluid system + FluidSystem::init(); + + // set the depth of the bottom of the reservoir + depthBOR_ = this->bBoxMax()[dimWorld-1]; + } + + //! Specifies the problem name. This is used as a prefix for files + //! generated by the simulation. + std::string name() const + { + static const std::string name = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, std::string, Problem, Name); + return name; + } + + //! Returns true if a restart file should be written. + bool shouldWriteRestartFile() const + { return false; } + + //! Returns true if the current solution should be written to disk + //! as a VTK file + bool shouldWriteOutput() const + { + return this->timeManager().timeStepIndex() > 0 && + this->timeManager().timeStepIndex() % 1 == 0; + } + + //! Returns the temperature within a finite volume. We use constant + //! 10 degrees Celsius. + Scalar temperature() const + { return 283.15; } + + //! Specifies which kind of boundary condition should be used for + //! which equation for a finite volume on the boundary. + void boundaryTypes(BoundaryTypes &bcTypes, const Vertex &vertex) const + { + const GlobalPosition &globalPos = vertex.geometry().center(); + + // Dirichlet conditions on left & right boundary + if (globalPos[0] < eps_ || globalPos[0] > this->bBoxMax()[0] - eps_) + bcTypes.setAllDirichlet(); + else // neuman for the remaining boundaries + bcTypes.setAllNeumann(); + + } + + //! Evaluates the Dirichlet boundary conditions for a finite volume + //! on the grid boundary. Here, the 'values' parameter stores + //! primary variables. + void dirichlet(PrimaryVariables &values, const Vertex &vertex) const + { + const auto globalPos = vertex.geometry().center(); + values[Indices::pressureIdx] = 200.0e3 + 9.81*1000*(depthBOR_ - globalPos[dimWorld-1]); + values[Indices::switchIdx] = 0.0; // 0 % oil saturation on left boundary + } + + //! Evaluates the boundary conditions for a Neumann boundary + //! segment. Here, the 'values' parameter stores the mass flux in + //! [kg/(m^2 * s)] in normal direction of each phase. Negative + //! values mean influx. + void neumann(PrimaryVariables &values, + const Element &element, + const FVElementGeometry &fvGeometry, + const Intersection &intersection, + int scvIdx, + int boundaryFaceIdx) const + { + const GlobalPosition &globalPos = + fvGeometry.boundaryFace[boundaryFaceIdx].ipGlobal; + Scalar up = this->bBoxMax()[dimWorld-1]; + // extraction of oil on the right boundary for approx. 1.e6 seconds + if (globalPos[dimWorld-1] > up - eps_ && globalPos[0] > 20 && globalPos[0] < 40) { + // oil outflux of 30 g/(m * s) on the right boundary. + // we solve for the mole balance, so we have to divide by the molar mass + values[Indices::contiWEqIdx] = 0; + values[Indices::contiNEqIdx] = -3e-2/FluidSystem::MyCompressibleComponent::molarMass(); + } else { + // no-flow on the remaining Neumann-boundaries. + values[Indices::contiWEqIdx] = 0; + values[Indices::contiNEqIdx] = 0; + } + } + + //! Evaluates the initial value for a control volume. For this + //! method, the 'values' parameter stores primary variables. + void initial(PrimaryVariables &values, + const Element &element, + const FVElementGeometry &fvGeometry, + int scvIdx) const + { + const auto globalPos = fvGeometry.subContVol[scvIdx].global; + values[Indices::pressureIdx] = 200.0e3 + 9.81*1000*(depthBOR_ - globalPos[dimWorld-1]); // 200 kPa = 2 bar + values[Indices::switchIdx] = 0.0; + } + + //! Evaluates the source term for all phases within a given + //! sub-control-volume. In this case, the 'values' parameter + //! stores the rate mass generated or annihilated per volume unit + //! in [kg / (m^3 * s)]. Positive values mean that mass is created. + void source(PrimaryVariables &values, + const Element &element, + const FVElementGeometry &fvGeometry, + int scvIdx) const + { + values[Indices::contiWEqIdx] = 0.0; + values[Indices::contiNEqIdx]= 0.0; + } + + //! Evaluates the initial phase presence + int initialPhasePresenceAtPos(const GlobalPosition &globalPos) const + { return Indices::wPhaseOnly; } + + /*! + * \brief Append all quantities of interest which can be derived + * from the solution of the current time step to the VTK + * writer. Adjust this in case of anisotropic permeabilities. + */ + void addOutputVtkFields() + { + // get the number of elements in the grid + unsigned numElements = this->gridView().size(/*codim=*/0); + + // create the scalar field required for the output of the lenses + typedef Dune::BlockVector > ScalarField; + ScalarField *isInLens = this->resultWriter().allocateManagedBuffer(numElements); + + // add data to the scalar field + for (const auto& element : elements(this->gridView())) + (*isInLens)[this->model().elementMapper().index(element)] = this->spatialParams().isInLens(element.geometry().center()); + + // attach to writer + this->resultWriter().attachCellData(*isInLens, "isInLens"); + } + +private: + // small epsilon value + Scalar eps_; + + // depth at the bottom of the reservoir + Scalar depthBOR_; +}; +} + +#endif diff --git a/tutorial/ex3/fluidsystems/h2omycompressiblecomponent.hh b/tutorial/ex3/fluidsystems/h2omycompressiblecomponent.hh new file mode 100644 index 0000000000000000000000000000000000000000..b02bce33982c161b01176f0ab74cf283891eb952 --- /dev/null +++ b/tutorial/ex3/fluidsystems/h2omycompressiblecomponent.hh @@ -0,0 +1,457 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/***************************************************************************** + * See the file COPYING for full copying permissions. * + * * + * This program is free software: you can redistribute it and/or modify * + * it under the terms of the GNU General Public License as published by * + * the Free Software Foundation, either version 2 of the License, or * + * (at your option) any later version. * + * * + * This program is distributed in the hope that it will be useful, * + * but WITHOUT ANY WARRANTY; without even the implied warranty of * + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * + * GNU General Public License for more details. * + * * + * You should have received a copy of the GNU General Public License * + * along with this program. If not, see . * + *****************************************************************************/ +/*! + * \file + * + * \brief A fluid system with water and a ficitious component, which is to be + * implemented in tutorial exercise 3a, as phases and components. This + * fluid system is to be implemented in exercise 3b. + */ +#ifndef DUMUX_H2O_MYCOMPRESSIBLECOMPONENT_FLUID_SYSTEM_HH +#define DUMUX_H2O_MYCOMPRESSIBLECOMPONENT_FLUID_SYSTEM_HH + +#include +#include + +#include + +// the ficitious component that was created in exercise 3a +#include + +// the binary coefficients corresponding to this fluid system +#include + +namespace Dumux +{ +namespace FluidSystems +{ + +/*! + * \brief A compositional fluid consisting of two liquid phases, + * which are water and a ficitious component from tutorial exercise 3a. + */ +template > > +class H2OMyCompressibleComponent + : public BaseFluidSystem< Scalar, H2OMyCompressibleComponent > +{ + typedef H2OMyCompressibleComponent ThisType; + typedef BaseFluidSystem Base; + +public: + typedef Dumux::MyCompressibleComponent MyCompressibleComponent; + typedef H2OType H2O; + + static const int numPhases = 2; + static const int numComponents = 2; + + static const int wPhaseIdx = 0; // index of the water phase + static const int nPhaseIdx = 1; // index of the NAPL phase + + static const int H2OIdx = 0; + static const int NAPLIdx = 1; + + // export component indices to indicate the main component + // of the corresponding phase at atmospheric pressure 1 bar + // and room temperature 20°C: + static const int wCompIdx = H2OIdx; + static const int nCompIdx = NAPLIdx; + + /*! + * \brief Initialize the fluid system's static parameters generically + * + * If a tabulated H2O component is used, we do our best to create + * tables that always work. + */ + static void init() + { + init(/*tempMin=*/273.15, + /*tempMax=*/623.15, + /*numTemp=*/100, + /*pMin=*/0.0, + /*pMax=*/20e6, + /*numP=*/200); + } + + /*! + * \brief Initialize the fluid system's static parameters using + * problem specific temperature and pressure ranges + * + * \param tempMin The minimum temperature used for tabulation of water [K] + * \param tempMax The maximum temperature used for tabulation of water [K] + * \param nTemp The number of ticks on the temperature axis of the table of water + * \param pressMin The minimum pressure used for tabulation of water [Pa] + * \param pressMax The maximum pressure used for tabulation of water [Pa] + * \param nPress The number of ticks on the pressure axis of the table of water + */ + static void init(Scalar tempMin, Scalar tempMax, unsigned nTemp, + Scalar pressMin, Scalar pressMax, unsigned nPress) + { + if (H2O::isTabulated) { + std::cout << "Initializing tables for the H2O fluid properties (" + << nTemp*nPress + << " entries).\n"; + + H2O::init(tempMin, tempMax, nTemp, + pressMin, pressMax, nPress); + } + } + + + /*! + * \brief Return whether a phase is liquid + * + * \param phaseIdx The index of the fluid phase to consider + */ + static bool isLiquid(int phaseIdx) + { + assert(0 <= phaseIdx && phaseIdx < numPhases); + return true; + } + + static bool isIdealGas(int phaseIdx) + { return H2O::gasIsIdeal() && MyCompressibleComponent::gasIsIdeal(); } + + /*! + * \brief Returns true if and only if a fluid phase is assumed to + * be an ideal mixture. + * + * We define an ideal mixture as a fluid phase where the fugacity + * coefficients of all components times the pressure of the phase + * are indepent on the fluid composition. This assumtion is true + * if Henry's law and Raoult's law apply. If you are unsure what + * this function should return, it is safe to return false. The + * only damage done will be (slightly) increased computation times + * in some cases. + * + * \param phaseIdx The index of the fluid phase to consider + */ + static bool isIdealMixture(int phaseIdx) + { + assert(0 <= phaseIdx && phaseIdx < numPhases); + return true; + } + + /*! + * \brief Returns true if and only if a fluid phase is assumed to + * be compressible. + * + * Compressible means that the partial derivative of the density + * to the fluid pressure is always larger than zero. + * + * \param phaseIdx The index of the fluid phase to consider + */ + static bool isCompressible(int phaseIdx) + { + assert(0 <= phaseIdx && phaseIdx < numPhases); + if (phaseIdx == wPhaseIdx) + // the water component decides for the water phase... + return H2O::liquidIsCompressible(); + + // the NAPL component decides for the napl phase... + return MyCompressibleComponent::liquidIsCompressible(); + } + + /*! + * \brief Return the human readable name of a phase (used in indices) + */ + static std::string phaseName(int phaseIdx) + { + switch (phaseIdx) { + case wPhaseIdx: return "w"; + case nPhaseIdx: return "n"; + }; + DUNE_THROW(Dune::InvalidStateException, "Invalid phase index " << phaseIdx); + } + + /*! + * \brief Return the human readable name of a component (used in indices) + */ + static std::string componentName(int compIdx) + { + switch (compIdx) { + case H2OIdx: return H2O::name(); + case NAPLIdx: return MyCompressibleComponent::name(); + }; + DUNE_THROW(Dune::InvalidStateException, "Invalid component index " << compIdx); + } + + /*! + * \brief Return the molar mass of a component in [kg/mol]. + */ + static Scalar molarMass(int compIdx) + { + switch (compIdx) { + case H2OIdx: return H2O::molarMass(); + case NAPLIdx: return MyCompressibleComponent::molarMass(); + }; + DUNE_THROW(Dune::InvalidStateException, "Invalid component index " << compIdx); + } + + /*! + * \brief Given all mole fractions in a phase, return the phase + * density [kg/m^3]. + */ + using Base::density; + template + static Scalar density(const FluidState &fluidState, int phaseIdx) + { + assert(0 <= phaseIdx && phaseIdx < numPhases); + if (phaseIdx == wPhaseIdx) { + // See: doctoral thesis of Steffen Ochs 2007 + // Steam injection into saturated porous media : process analysis including experimental and numerical investigations + // http://elib.uni-stuttgart.de/bitstream/11682/271/1/Diss_Ochs_OPUS.pdf + Scalar rholH2O = H2O::liquidDensity(fluidState.temperature(phaseIdx), fluidState.pressure(phaseIdx)); + Scalar clH2O = rholH2O/H2O::molarMass(); + Scalar x_H2O = fluidState.moleFraction(wPhaseIdx, H2OIdx); + Scalar x_myComp = fluidState.moleFraction(wPhaseIdx, NAPLIdx); + + /*! + * TODO: implement the composition-dependent water density from the exercise sheet. + */ + return ???; + } + else { + // assume the density of the fictious component to be independent of the composition + Scalar pressure = MyCompressibleComponent::liquidIsCompressible()?fluidState.pressure(phaseIdx):1e100; + return MyCompressibleComponent::liquidDensity(fluidState.temperature(phaseIdx), pressure); + } + } + + /*! + * \brief Return the viscosity of a phase. + */ + using Base::viscosity; + template + static Scalar viscosity(const FluidState &fluidState, + int phaseIdx) + { + assert(0 <= phaseIdx && phaseIdx < numPhases); + if (phaseIdx == wPhaseIdx) { + // assume pure water viscosity + return H2O::liquidViscosity(fluidState.temperature(phaseIdx), + fluidState.pressure(phaseIdx)); + } + else { + // assume pure NAPL viscosity + return MyCompressibleComponent::liquidViscosity(fluidState.temperature(phaseIdx), fluidState.pressure(phaseIdx)); + } + } + + using Base::diffusionCoefficient; + template + static Scalar diffusionCoefficient(const FluidState &fluidState, int phaseIdx, int compIdx) + { + DUNE_THROW(Dune::NotImplemented, "Diffusion coefficients"); + } + + /*! + * \brief Given a phase's composition, temperature and pressure, + * return the binary diffusion coefficient \f$\mathrm{[m^2/s]}\f$ for components + * \f$\mathrm{i}\f$ and \f$\mathrm{j}\f$ in this phase. + * \param fluidState The fluid state + * \param paramCache mutable parameters + * \param phaseIdx Index of the fluid phase + * \param compIIdx Index of the component i + * \param compJIdx Index of the component j + */ + using Base::binaryDiffusionCoefficient; + template + static Scalar binaryDiffusionCoefficient(const FluidState &fluidState, + int phaseIdx, + int compIIdx, + int compJIdx) + { + assert(0 <= phaseIdx && phaseIdx < numPhases); + assert(0 <= compIdx && compIdx < numComponents); + + const Scalar T = fluidState.temperature(phaseIdx); + const Scalar p = fluidState.pressure(phaseIdx); + + // we assume the diffusion coefficient to be the same in both phases + return Dumux::BinaryCoeff::H2O_MyCompressibleComponent::liquidDiffCoeff(T, p); + } + + /* Henry coefficients + */ + template + static Scalar henryCoefficient(const FluidState &fluidState, + int phaseIdx, + int compIdx) + { + assert(0 <= phaseIdx && phaseIdx < numPhases); + assert(0 <= compIdx && compIdx < numComponents); + + const Scalar T = fluidState.temperature(phaseIdx); + const Scalar p = fluidState.pressure(phaseIdx); + + if (compIdx == NAPLIdx && phaseIdx == wPhaseIdx) + return Dumux::BinaryCoeff::H2O_MyCompressibleComponent::henryMyCompressibleComponentInWater(T)/p; + + else if (phaseIdx == nPhaseIdx && compIdx == H2OIdx) + return Dumux::BinaryCoeff::H2O_MyCompressibleComponent::henryWaterInMyCompressibleComponent(T)/p; + + else + DUNE_THROW(Dune::InvalidStateException, "non-existent henry coefficient for phase index " << phaseIdx + << " and component index " << compIdx); + } + + using Base::fugacityCoefficient; + /*! + * \brief Returns the fugacity coefficient \f$\mathrm{[-]}\f$ of a component in a + * phase. + * + * In this case, things are actually pretty simple. We have an ideal + * solution. Thus, the fugacity coefficient is 1 in the gas phase + * (fugacity equals the partial pressure of the component in the gas phase + * respectively in the liquid phases it is the inverse of the + * Henry coefficients scaled by pressure + * \param fluidState The fluid state + * \param phaseIdx The index of the phase + * \param compIdx The index of the component + */ + template + static Scalar fugacityCoefficient(const FluidState &fluidState, + int phaseIdx, + int compIdx) + { + assert(0 <= phaseIdx && phaseIdx < numPhases); + assert(0 <= compIdx && compIdx < numComponents); + + Scalar T = fluidState.temperature(phaseIdx); + Scalar p = fluidState.pressure(phaseIdx); + + if (phaseIdx == wPhaseIdx) { + if (compIdx == H2OIdx) + return H2O::vaporPressure(T)/p; + else if (compIdx == NAPLIdx) + return Dumux::BinaryCoeff::H2O_MyCompressibleComponent::henryMyCompressibleComponentInWater(T)/p; + } + + // for the NAPL phase, we assume currently that nothing is + // dissolved. this means that the affinity of the NAPL + // component to the NAPL phase is much higher than for the + // other components, i.e. the fugacity coefficient is much + // smaller. + Scalar phiNapl = MyCompressibleComponent::vaporPressure(T)/p; + if (compIdx == NAPLIdx) + return phiNapl; + else + return 1e6*phiNapl; + } + + template + static Scalar kelvinVaporPressure(const FluidState &fluidState, + const int phaseIdx, + const int compIdx) + { + DUNE_THROW(Dune::NotImplemented, "FluidSystems::H2OMyCompressibleComponent::kelvinVaporPressure()"); + } + + /* partial pressures in the gas phase, taken from saturation vapor pressures + */ + template + static Scalar partialPressureGas(const FluidState &fluidState, int phaseIdx, + int compIdx) + { + assert(0 <= compIdx && compIdx < numComponents); + + const Scalar T = fluidState.temperature(phaseIdx); + if (compIdx == NAPLIdx) + return MyCompressibleComponent::vaporPressure(T); + else if (compIdx == H2OIdx) + return H2O::vaporPressure(T); + else + DUNE_THROW(Dune::InvalidStateException, "non-existent component index " << compIdx); + } + + /* inverse vapor pressures, taken from inverse saturation vapor pressures + */ + template + static Scalar inverseVaporPressureCurve(const FluidState &fluidState, + int phaseIdx, + int compIdx) + { + assert(0 <= compIdx && compIdx < numComponents); + + const Scalar pressure = fluidState.pressure(phaseIdx); + if (compIdx == NAPLIdx) + return MyCompressibleComponent::vaporTemperature(pressure); + else if (compIdx == H2OIdx) + return H2O::vaporTemperature(pressure); + else + DUNE_THROW(Dune::InvalidStateException, "non-existent component index " << compIdx); + } + + + + /*! + * \brief Given all mole fractions in a phase, return the specific + * phase enthalpy [J/kg]. + */ + using Base::enthalpy; + template + static Scalar enthalpy(const FluidState &fluidState, + int phaseIdx) + { + assert(0 <= phaseIdx && phaseIdx < numPhases); + if (phaseIdx == wPhaseIdx) { + return H2O::liquidEnthalpy(fluidState.temperature(phaseIdx), fluidState.pressure(phaseIdx)); + } + else { + return MyCompressibleComponent::liquidEnthalpy(fluidState.temperature(phaseIdx), fluidState.pressure(phaseIdx)); + } + DUNE_THROW(Dune::InvalidStateException, "Invalid phase index " << phaseIdx); + } + + using Base::heatCapacity; + template + static Scalar heatCapacity(const FluidState &fluidState, + int phaseIdx) + { + DUNE_THROW(Dune::NotImplemented, "FluidSystems::H2ONAPL::heatCapacity()"); + } + + using Base::thermalConductivity; + template + static Scalar thermalConductivity(const FluidState &fluidState, + int phaseIdx) + { + assert(0 <= phaseIdx && phaseIdx < numPhases); + + const Scalar temperature = fluidState.temperature(phaseIdx) ; + const Scalar pressure = fluidState.pressure(phaseIdx); + if (phaseIdx == wPhaseIdx) + { + return H2O::liquidThermalConductivity(temperature, pressure); + } + else + { + return MyCompressibleComponent::liquidThermalConductivity(temperature, pressure); + } + DUNE_THROW(Dune::InvalidStateException, "Invalid phase index " << phaseIdx); + } + +private: + +}; +} // end namespace FluidSystems +} // end namespace Dumux + +#endif diff --git a/tutorial/ex3/spatialparams.hh b/tutorial/ex3/spatialparams.hh new file mode 100644 index 0000000000000000000000000000000000000000..4d5f0635f99d8ca0a426f67d0f1c42fdb3102f50 --- /dev/null +++ b/tutorial/ex3/spatialparams.hh @@ -0,0 +1,197 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/***************************************************************************** + * See the file COPYING for full copying permissions. * + * * + * This program is free software: you can redistribute it and/or modify * + * it under the terms of the GNU General Public License as published by * + * the Free Software Foundation, either version 2 of the License, or * + * (at your option) any later version. * + * * + * This program is distributed in the hope that it will be useful, * + * but WITHOUT ANY WARRANTY; without even the implied warranty of * + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * + * GNU General Public License for more details. * + * * + * You should have received a copy of the GNU General Public License * + * along with this program. If not, see . * + *****************************************************************************/ +/*! + * \file + * + * \brief The spatial parameters for the fully coupled tutorial problem + * which uses the twophase box model. + */ +#ifndef DUMUX_EXERCISE_THREE_SPATIAL_PARAMS_HH +#define DUMUX_EXERCISE_THREE_SPATIAL_PARAMS_HH + +// include parent spatialparameters +#include + +// include material laws +#include +#include +#include + +namespace Dumux { +//forward declaration +template +class ExerciseThreeSpatialParams; + +namespace Properties +{ +// The spatial parameters TypeTag +NEW_TYPE_TAG(ExerciseThreeSpatialParams); + +// Set the spatial parameters +SET_TYPE_PROP(ExerciseThreeSpatialParams, SpatialParams, + ExerciseThreeSpatialParams); + +// Set the material law +SET_PROP(ExerciseThreeSpatialParams, MaterialLaw) +{ +private: + // material law typedefs + typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar; + // select material law to be used + typedef RegularizedBrooksCorey RawMaterialLaw; +public: + // adapter for absolute law + typedef EffToAbsLaw type; +}; +} + +/*! + * \ingroup TwoPBoxModel + * + * \brief The spatial parameters for the fully coupled tutorial problem + * which uses the twophase box model. + */ +template +class ExerciseThreeSpatialParams: public ImplicitSpatialParams +{ + // Get informations for current implementation via property system + typedef typename GET_PROP_TYPE(TypeTag, Grid) Grid; + typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView; + typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar; + enum + { + dim = Grid::dimension, + dimWorld = Grid::dimensionworld + }; + + // Get object types for function arguments + typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry; + typedef typename Grid::Traits::template Codim<0>::Entity Element; + +public: + // get material law from property system + typedef typename GET_PROP_TYPE(TypeTag, MaterialLaw) MaterialLaw; + // determine appropriate parameters depending on selected materialLaw + typedef typename MaterialLaw::Params MaterialLawParams; + + /*! Intrinsic permeability tensor K \f$[m^2]\f$ depending + * on the position in the domain + * + * \param element The finite volume element + * \param fvGeometry The finite-volume geometry in the box scheme + * \param scvIdx The local vertex index + * + * Alternatively, the function intrinsicPermeabilityAtPos(const GlobalPosition& globalPos) + * could be defined, where globalPos is the vector including the global coordinates + * of the finite volume. + */ + const Dune::FieldMatrix &intrinsicPermeability(const Element &element, + const FVElementGeometry &fvGeometry, + const int scvIdx) const + { + if (isInLens(element.geometry().center())) + return KLens_; + return K_; + } + + /*! Defines the porosity \f$[-]\f$ of the porous medium depending + * on the position in the domain + * + * \param element The finite volume element + * \param fvGeometry The finite-volume geometry in the box scheme + * \param scvIdx The local vertex index + * + * Alternatively, the function porosityAtPos(const GlobalPosition& globalPos) + * could be defined, where globalPos is the vector including the global coordinates + * of the finite volume. + */ + Scalar porosity(const Element &element, + const FVElementGeometry &fvGeometry, + const int scvIdx) const + { + if (isInLens(element.geometry().center())) + return 0.1; + return 0.2; + } + + /*! Returns the parameter object for the material law (i.e. Brooks-Corey) + * depending on the position in the domain + * + * \param element The finite volume element + * \param fvGeometry The finite-volume geometry in the box scheme + * \param scvIdx The local vertex index + * + * Alternatively, the function materialLawParamsAtPos(const GlobalPosition& globalPos) + * could be defined, where globalPos is the vector including the global coordinates + * of the finite volume. + */ + const MaterialLawParams& materialLawParams(const Element &element, + const FVElementGeometry &fvGeometry, + const int scvIdx) const + { + if (isInLens(element.geometry().center())) + return materialParamsLens_; + return materialParams_; + } + + // constructor + ExerciseThreeSpatialParams(const GridView& gridView) : + ImplicitSpatialParams(gridView), + K_(0), + KLens_(0) + { + //set main diagonal entries of the permeability tensor to a value + //setting to one value means: isotropic, homogeneous + for (int i = 0; i < dim; i++) + { + K_[i][i] = 1e-7; + KLens_[i][i] = 1e-10; + } + + //set residual saturations + materialParams_.setSwr(0.0); + materialParamsLens_.setSwr(0.1); + materialParams_.setSnr(0.0); + materialParamsLens_.setSnr(0.1); + + //parameters of Brooks & Corey Law + materialParams_.setPe(500.0); + materialParamsLens_.setPe(1000.0); + materialParams_.setLambda(2); + materialParamsLens_.setLambda(2); + } + + bool isInLens(const Dune::FieldVector& globalPos) const + { + const auto x = globalPos[0]; + const auto y = globalPos[1]; + return (x < 40 && x > 20 && y > 35 && y < 45) || + (x < 50 && x > 30 && y < 30 && y > 15); + } + +private: + + Dune::FieldMatrix K_; + Dune::FieldMatrix KLens_; + // Object that holds the values/parameters of the selected material law. + MaterialLawParams materialParams_; + MaterialLawParams materialParamsLens_; +}; +} // end namespace +#endif diff --git a/tutorial/ex4/README.md b/tutorial/ex4/README.md new file mode 100644 index 0000000000000000000000000000000000000000..de195573c2340ed3719405fb14e877ff164ce1e8 --- /dev/null +++ b/tutorial/ex4/README.md @@ -0,0 +1,92 @@ +# Exercise #4 (DuMuX course) + +This exercise describes how to create a new DuMuX module +and how to create a corresponding GitLab project. + +This is the suggested +workflow to develop code on top of DuMuX. + +### Task 1: Create new dune module +
+ +* Execute the following command (bash environment) in the top-folder, i.e. above the dumux folder + +```bash +./dune-common/bin/duneproject +``` + +* Follow the introductions and specify + * as name of the new module: `dumux-example` + * as module dependencies: `dumux` + * a version at your choice + * your email adress + +


+### Task 2: Rerun dunecontrol to configure your new project +
+ +The following command will configure your new module + +```bash +./dune-common/bin/dunecontrol --opts= --only=dumux-example all +``` + +


+### Task 3: Create a new test case within your new DuMuX module +
+ +* Create a new folder (in your module folder), e.g. `appl` + +```bash +mkdir appl +``` + +* Copy some test case from the dumux module, e.g. test_box1p + * Copy the problem, spatialparams, cc source file, input file + +* Adjust the CMakeLists.txt file to include your new subdirectory + +* Add a new CMakeLists.txt in the folder `appl` with the content + +```cmake +// add a new box 1p test +dune_add_test(NAME test_box1p + SOURCES test_box1p.cc) + +// link the input file to the build folder +dune_symlink_to_source_files(FILES test_box1p.input) + +``` + +* Reconfigure your module by running in the topmost directory of your new module + +```bash +cmake build-cmake +``` + +* Build and execute the test problem + +```bash +cd build-cmake +make build_tests +cd appl +./test_box1p +``` + +


+### Task 4: Create a new gitlab project +
+ +* Login with your username and password at https://git.iws.uni-stuttgart.de/ + +Note: If you don't have an account create one. We allow anyone to host repositories +on our GitLab instance as long as it is DuMuX related. + +* Click the **New project** button + +* Follow the given instructions for an *existing folder* + +**Important**: Before executing the `git add .` command, you should add your cmake build folder to `.gitignore`. +The easiest way to do so is to copy the `.gitignore` file from the dumux module into your module path. If everything +worked, executing `git status` should not show `build-cmake` anymore. Never put your executables or other build files +under version control. Only source files (`*.hh`, `*.cc`, `*.input`, `CMakeLists.txt`) should be under version control. diff --git a/tutorial/extradoc/exercise1_nonisothermal.png b/tutorial/extradoc/exercise1_nonisothermal.png new file mode 100644 index 0000000000000000000000000000000000000000..938ba158eb97e83d1ae80145161b94d2ab6137cf Binary files /dev/null and b/tutorial/extradoc/exercise1_nonisothermal.png differ diff --git a/tutorial/extradoc/exercise1_setup.png b/tutorial/extradoc/exercise1_setup.png new file mode 100644 index 0000000000000000000000000000000000000000..db9c1b4acdfe9368781cd7593b0def3963146844 Binary files /dev/null and b/tutorial/extradoc/exercise1_setup.png differ diff --git a/tutorial/extradoc/exercise2_properties.png b/tutorial/extradoc/exercise2_properties.png new file mode 100644 index 0000000000000000000000000000000000000000..28c700ed18443d78995d890cb7fb085a84146ff7 Binary files /dev/null and b/tutorial/extradoc/exercise2_properties.png differ diff --git a/tutorial/extradoc/exercise3_a_solution.png b/tutorial/extradoc/exercise3_a_solution.png new file mode 100644 index 0000000000000000000000000000000000000000..fdf3ef23c81a0bc414554a9bc2e4f1e784e134c4 Binary files /dev/null and b/tutorial/extradoc/exercise3_a_solution.png differ diff --git a/tutorial/extradoc/exercise3_a_solution2.png b/tutorial/extradoc/exercise3_a_solution2.png new file mode 100644 index 0000000000000000000000000000000000000000..ba90e6fcfa783fcc5a4eaf1c5fbe252dc46ea8d9 Binary files /dev/null and b/tutorial/extradoc/exercise3_a_solution2.png differ diff --git a/tutorial/extradoc/exercise3_setup.png b/tutorial/extradoc/exercise3_setup.png new file mode 100644 index 0000000000000000000000000000000000000000..eed8371e6a96f9b4d4eee22849f4b3bbbcd83fa1 Binary files /dev/null and b/tutorial/extradoc/exercise3_setup.png differ diff --git a/tutorial/solution/ex1/exercise1_2pni.cc b/tutorial/solution/ex1/exercise1_2pni.cc new file mode 100644 index 0000000000000000000000000000000000000000..69f233e059a5c2aba48928a415df905dec239e45 --- /dev/null +++ b/tutorial/solution/ex1/exercise1_2pni.cc @@ -0,0 +1,64 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/***************************************************************************** + * See the file COPYING for full copying permissions. * + * * + * This program is free software: you can redistribute it and/or modify * + * it under the terms of the GNU General Public License as published by * + * the Free Software Foundation, either version 2 of the License, or * + * (at your option) any later version. * + * * + * This program is distributed in the hope that it will be useful, * + * but WITHOUT ANY WARRANTY; without even the implied warranty of * + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * + * GNU General Public License for more details. * + * * + * You should have received a copy of the GNU General Public License * + * along with this program. If not, see . * + *****************************************************************************/ +/*! + * \file + * + * \brief DOC ME! + */ +#include "config.h" +#include "injection2pniproblem.hh" +#include + +/*! + * \brief Provides an interface for customizing error messages associated with + * reading in parameters. + * + * \param progName The name of the program, that was tried to be started. + * \param errorMsg The error message that was issued by the start function. + * Comprises the thing that went wrong and a general help message. + */ +void usage(const char *progName, const std::string &errorMsg) +{ + if (errorMsg.size() > 0) { + std::string errorMessageOut = "\nUsage: "; + errorMessageOut += progName; + errorMessageOut += " [options]\n"; + errorMessageOut += errorMsg; + errorMessageOut += "\n\nThe list of mandatory options for this program is:\n" + "\t-TimeManager.TEnd End of the simulation [s] \n" + "\t-TimeManager.DtInitial Initial timestep size [s] \n" + "\t-Grid.File Name of the file containing the grid \n" + "\t definition in DGF format\n" + "\t-SpatialParams.LensLowerLeft coordinates of the lower left corner of the lens [m] \n" + "\t-SpatialParams.LensUpperRight coordinates of the upper right corner of the lens [m] \n" + "\t-Problem.Name String for naming of the output files \n" + "\n"; + std::cout << errorMessageOut + << "\n"; + } +} + +//////////////////////// +// the main function +//////////////////////// +int main(int argc, char** argv) +{ + typedef TTAG(InjectionCCProblem2PNI) ProblemTypeTag; + return Dumux::start(argc, argv, usage); +} diff --git a/tutorial/solution/ex1/injection2p2cproblem.hh b/tutorial/solution/ex1/injection2p2cproblem.hh new file mode 100644 index 0000000000000000000000000000000000000000..e03f339f15072819f16f91b52a31af1f75173640 --- /dev/null +++ b/tutorial/solution/ex1/injection2p2cproblem.hh @@ -0,0 +1,310 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/***************************************************************************** + * See the file COPYING for full copying permissions. * + * * + * This program is free software: you can redistribute it and/or modify * + * it under the terms of the GNU General Public License as published by * + * the Free Software Foundation, either version 2 of the License, or * + * (at your option) any later version. * + * * + * This program is distributed in the hope that it will be useful, * + * but WITHOUT ANY WARRANTY; without even the implied warranty of * + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * + * GNU General Public License for more details. * + * * + * You should have received a copy of the GNU General Public License * + * along with this program. If not, see . * + *****************************************************************************/ +/*! + * \file + * + * \brief Problem where air is injected under a low permeable layer in a depth of 2700m. + */ +#ifndef DUMUX_INJECTION_2P2C_PROBLEM_HH +#define DUMUX_INJECTION_2P2C_PROBLEM_HH + +#include +#include +#include + +#include "injection2p2cspatialparams.hh" + +namespace Dumux +{ + +// foward declaration +template +class Injection2p2cProblem; + +// setup property TypeTag +namespace Properties +{ + +// inherit from MyLocalResidualParams +NEW_TYPE_TAG(Injection2p2cProblem, INHERITS_FROM(TwoPTwoC, InjectionSpatialParams)); +NEW_TYPE_TAG(Injection2p2cBoxProblem, INHERITS_FROM(BoxModel, Injection2p2cProblem)); +NEW_TYPE_TAG(Injection2p2pcCCProblem, INHERITS_FROM(CCModel, Injection2p2cProblem)); + +// Set the grid type +SET_TYPE_PROP(Injection2p2cProblem, Grid, Dune::YaspGrid<2>); + +// Set the problem property +SET_TYPE_PROP(Injection2p2cProblem, Problem, Injection2p2cProblem); + +// Set fluid configuration +SET_TYPE_PROP(Injection2p2cProblem, + FluidSystem, + FluidSystems::H2ON2); + +// Define whether mole(true) or mass (false) fractions are used +SET_BOOL_PROP(Injection2p2cProblem, UseMoles, true); + +} // end namespace Properties + +/*! + * \ingroup TwoPTwoCModel + * \ingroup ImplicitTestProblems + * \brief Problem where air is injected under a low permeable layer in a depth of 2700m. + * + * The domain is sized 60m times 40m and consists of two layers, a moderately + * permeable one for \f$ y<22m\f$ and one with a lower permeablility + * in the rest of the domain. + * + * Nitrogen is injected into a water-filled aquifer through a well. First, we inject for one month. + * Then, we continue simulating the development of the nitrogen plume for 10 years. + * This is realized with a Neumann boundary condition at the right boundary + * (\f$ 7m +class Injection2p2cProblem : public ImplicitPorousMediaProblem +{ + typedef ImplicitPorousMediaProblem ParentType; + + typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView; + typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar; + typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem; + + // grid world dimension + static constexpr auto dimWorld = GridView::dimensionworld; + + typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices; + typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables; + typedef typename GET_PROP_TYPE(TypeTag, BoundaryTypes) BoundaryTypes; + typedef typename GET_PROP_TYPE(TypeTag, TimeManager) TimeManager; + typedef typename GET_PROP_TYPE(TypeTag, MaterialLaw) MaterialLaw; + typedef Dune::FieldVector GlobalPosition; + +public: + /*! + * \brief The constructor + * + * \param timeManager The time manager + * \param gridView The grid view + */ + Injection2p2cProblem(TimeManager &timeManager, const GridView &gridView) + : ParentType(timeManager, gridView) + { + // initialize the tables of the fluid system + FluidSystem::init(/*tempMin=*/273.15, + /*tempMax=*/423.15, + /*numTemp=*/50, + /*pMin=*/0.0, + /*pMax=*/30e6, + /*numP=*/300); + + // name of the problem and output file + name_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, std::string, Problem, Name); + // depth of the aquifer, units: m + aquiferDepth_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Problem, AquiferDepth); + // inflow rate of nitrogen water vapor mixture, units: kg/(s m^2) + totalAreaSpecificInflow_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Problem, TotalAreaSpecificInflow); + // the duration of the injection, units: second + injectionDuration_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Problem, InjectionDuration); + } + + /*! + * \brief User defined output after the time integration + * + * Will be called diretly after the time integration + */ + void postTimeStep() + { + // Calculate storage terms + PrimaryVariables storageW, storageN; + this->model().globalPhaseStorage(storageW, Indices::wPhaseIdx); + this->model().globalPhaseStorage(storageN, Indices::nPhaseIdx); + + // Write mass balance information for rank 0 + if (this->gridView().comm().rank() == 0) + { + std::cout <<"Storage: wetting=[" << storageW << "]" + << " nonwetting=[" << storageN << "]" << std::endl; + } + } + + /*! + * \name Problem parameters + */ + // \{ + + /*! + * \brief Returns the problem name + * + * This is used as a prefix for files generated by the simulation. + */ + const std::string& name() const + { return name_; } + + /*! + * \brief Returns the temperature in \f$ K \f$ + */ + Scalar temperature() const + { return 273.15 + 30; } + + /*! + * \brief Returns the source term + * + * \param values Stores the source values for the conservation equations in + * \f$ [ \textnormal{unit of primary variable} / (m^\textrm{dim} \cdot s )] \f$ + * \param globalPos The global position + */ + void sourceAtPos(PrimaryVariables &values, + const GlobalPosition &globalPos) const + { values = 0; } + + // \} + + /*! + * \name Boundary conditions + */ + // \{ + + /*! + * \brief Specifies which kind of boundary condition should be + * used for which equation on a given boundary segment + * + * \param values Stores the value of the boundary type + * \param globalPos The global position + */ + void boundaryTypesAtPos(BoundaryTypes &values, + const GlobalPosition &globalPos) const + { + // Set Dirichlet at the bottom of the domain + if (globalPos[dimWorld-1] < eps_) + { + values.setAllDirichlet(); + } + + // and Neuman boundary conditions everywhere else + // note that we don't differentiate between Neumann and Robin boundary types + else + { + values.setAllNeumann(); + } + } + + /*! + * \brief Evaluates the boundary conditions for a Dirichlet + * boundary segment + * + * \param values Stores the Dirichlet values for the conservation equations in + * \f$ [ \textnormal{unit of primary variable} ] \f$ + * \param globalPos The global position + */ + void dirichletAtPos(PrimaryVariables &values, const GlobalPosition &globalPos) const + { initialAtPos(values, globalPos); } + + /*! + * \brief Evaluates the boundary conditions for a Neumann boundary segment. + * + * \param values Stores the Neumann values for the conservation equations in + * \f$ [ \textnormal{unit of conserved quantity} / (m^(dim-1) \cdot s )] \f$ + * \param globalPos The globalPosition of the boundary interface + */ + void neumannAtPos(PrimaryVariables &values, + const GlobalPosition &globalPos) const + { + // initialize values to zero, i.e. no-flow Neumann boundary conditions + values = 0; + + //if we are inside the injection zone set inflow Neumann boundary conditions + if (this->timeManager().time() + this->timeManager().timeStepSize() < injectionDuration_ + && globalPos[1] < 15 + eps_ && globalPos[1] > 7 - eps_ && globalPos[0] > 0.9*this->bBoxMax()[0]) + { + // set the Neumann values for the Nitrogen component balance + // convert from units kg/(s*m^2) to mole/(s*m^2) + values[Indices::contiNEqIdx] = -totalAreaSpecificInflow_/FluidSystem::molarMass(FluidSystem::nCompIdx); + values[Indices::contiWEqIdx] = 0.0; + } + } + + // \} + + /*! + * \name Volume terms + */ + // \{ + + /*! + * \brief Evaluates the initial values for a control volume + * + * \param values Stores the initial values for the conservation equations in + * \f$ [ \textnormal{unit of primary variables} ] \f$ + * \param globalPos The global position + */ + void initialAtPos(PrimaryVariables &values, const GlobalPosition &globalPos) const + { + // get the water density at atmospheric conditions + const Scalar densityW = FluidSystem::H2O::liquidDensity(temperature(), 1.0e5); + + // assume an intially hydrostatic liquid pressure profile + // note: we subtract rho_w*g*h because g is defined negative + const Scalar pw = 1.0e5 - densityW*this->gravity()[dimWorld-1]*(aquiferDepth_ - globalPos[dimWorld-1]); + + // initially we have some nitrogen dissolved + // saturation mole fraction would be + // moleFracLiquidN2 = (pw + pc + p_vap^sat)/henry; + const Scalar moleFracLiquidN2 = pw*0.95/BinaryCoeff::H2O_N2::henry(temperature()); + + // note that because we start with a single phase system the primary variables + // are pl and x^w_N2. This will switch as soon after we start injecting to a two + // phase system so the primary variables will be pl and Sn (non-wetting saturation). + values[Indices::switchIdx] = moleFracLiquidN2; + values[Indices::pressureIdx] = pw; + } + + /*! + * \brief Return the initial phase state inside a control volume. + * + * \param globalPos The global position + * \note we start with a single phase system + */ + int initialPhasePresenceAtPos(const GlobalPosition &globalPos) const + { return Indices::wPhaseOnly; } + + // \} + + //! If we should write restart files + bool shouldWriteRestartFile() const + { return false; } + +private: + static constexpr Scalar eps_ = 1e-6; + std::string name_; //! Problem name + Scalar aquiferDepth_; //! Depth of the aquifer in m + Scalar totalAreaSpecificInflow_; //! Area specific inflow rate in mole/(s*m^2) + Scalar injectionDuration_; //! Duration of the injection in seconds +}; + +} // end namespace Dumux + +#endif diff --git a/tutorial/solution/ex1/injection2pniproblem.hh b/tutorial/solution/ex1/injection2pniproblem.hh new file mode 100644 index 0000000000000000000000000000000000000000..a023716ea2ea6408e6b627722d4148eb12352f05 --- /dev/null +++ b/tutorial/solution/ex1/injection2pniproblem.hh @@ -0,0 +1,291 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/***************************************************************************** + * See the file COPYING for full copying permissions. * + * * + * This program is free software: you can redistribute it and/or modify * + * it under the terms of the GNU General Public License as published by * + * the Free Software Foundation, either version 2 of the License, or * + * (at your option) any later version. * + * * + * This program is distributed in the hope that it will be useful, * + * but WITHOUT ANY WARRANTY; without even the implied warranty of * + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * + * GNU General Public License for more details. * + * * + * You should have received a copy of the GNU General Public License * + * along with this program. If not, see . * + *****************************************************************************/ +/*! + * \file + * + * \brief Non-isothermal gas injection problem where a gas (e.g. air) is injected into a fully + * water saturated medium. During buoyancy driven upward migration the gas + * passes a high temperature area. + */ + +#ifndef DUMUX_INJECTION_PROBLEM_2PNI_HH +#define DUMUX_INJECTION_PROBLEM_2PNI_HH + +#include +#include + +#include + +#include "injection2pspatialparams.hh" + +namespace Dumux { + +template +class InjectionProblem2PNI; + +namespace Properties +{ +NEW_TYPE_TAG(InjectionProblem2PNI, INHERITS_FROM(TwoPNI, InjectionSpatialParams)); +NEW_TYPE_TAG(InjectionCCProblem2PNI, INHERITS_FROM(CCModel, InjectionProblem2PNI)); + +// Set the grid type +SET_TYPE_PROP(InjectionProblem2PNI, Grid, Dune::YaspGrid<2>); + +// Set the problem property +SET_TYPE_PROP(InjectionProblem2PNI, Problem, InjectionProblem2PNI); + +// Use the same fluid system as the 2p2c injection problem +SET_TYPE_PROP(InjectionProblem2PNI, FluidSystem, FluidSystems::H2ON2); + +} + +/*! + * \ingroup TwoPModel + * \ingroup ImplicitTestProblems + * \brief Non-isothermal gas injection problem where a gas (e.g. air) is injected into a fully + * water saturated medium. During buoyancy driven upward migration the gas + * passes a high temperature area. + * + * The domain is sized 60 m times 40 m. The rectangular area with the increased temperature (380 K) + * starts at (20 m, 5 m) and ends at (30 m, 35 m) + * + * For the mass conservation equation neumann boundary conditions are used on + * the top, on the bottom and on the right of the domain, while dirichlet conditions + * apply on the left boundary. + * For the energy conservation equation dirichlet boundary conditions are applied + * on all boundaries. + * + * Gas is injected at the right boundary from 7 m to 15 m at a rate of + * 0.001 kg/(s m), the remaining neumann boundaries are no-flow + * boundaries. + * + * At the dirichlet boundaries a hydrostatic pressure, a gas saturation of zero and + * a geothermal temperature gradient of 0.03 K/m are applied. + * + * This problem uses the \ref TwoPModel and \ref NIModel model. + * + * To run the simulation execute the following line in shell: + * ./exercise1_2pni -ParameterFile exercise1.input + */ +template +class InjectionProblem2PNI : public ImplicitPorousMediaProblem +{ + typedef ImplicitPorousMediaProblem ParentType; + typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView; + typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar; + + + typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices; + enum { + pressureIdx = Indices::pressureIdx, + saturationIdx = Indices::saturationIdx, + + contiNEqIdx = Indices::contiNEqIdx, + + temperatureIdx = Indices::temperatureIdx, + energyEqIdx = Indices::energyEqIdx, + + // world dimension + dimWorld = GridView::dimensionworld + }; + + + typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables; + typedef typename GET_PROP_TYPE(TypeTag, BoundaryTypes) BoundaryTypes; + typedef typename GET_PROP_TYPE(TypeTag, TimeManager) TimeManager; + + typedef typename GridView::template Codim<0>::Entity Element; + typedef typename GridView::Intersection Intersection; + + typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry; + typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem; + + typedef Dune::FieldVector GlobalPosition; + + enum { isBox = GET_PROP_VALUE(TypeTag, ImplicitIsBox) }; + +public: + /*! + * \brief The constructor + * + * \param timeManager The time manager + * \param gridView The grid view + */ + InjectionProblem2PNI(TimeManager &timeManager, const GridView &gridView) + : ParentType(timeManager, gridView) + { + // initialize the tables of the fluid system + FluidSystem::init(/*tempMin=*/273.15, + /*tempMax=*/423.15, + /*numTemp=*/50, + /*pMin=*/0.0, + /*pMax=*/30e6, + /*numP=*/300); + + // name of the problem and output file + name_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, std::string, Problem, Name); + // depth of the aquifer, units: m + aquiferDepth_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Problem, AquiferDepth); + // inflow rate of nitrogen water vapor mixture, units: kg/(s m^2) + // the duration of the injection, units: second + injectionDuration_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Problem, InjectionDuration); + } + + /*! + * \name Problem parameters + */ + // \{ + + /*! + * \brief Returns the problem name + * + * This is used as a prefix for files generated by the simulation. + */ + const std::string name() const + { return "injection-2pni";} + + /*! + * \brief Returns the source term + * + * \param values Stores the source values for the conservation equations in + * \f$ [ \textnormal{unit of primary variable} / (m^\textrm{dim} \cdot s )] \f$ + * \param globalPos The global position + */ + void sourceAtPos(PrimaryVariables &values, + const GlobalPosition &globalPos) const + { + values = 0; + } + + // \} + + /*! + * \name Boundary conditions + */ + // \{ + + /*! + * \brief Specifies which kind of boundary condition should be + * used for which equation on a given boundary segment + * + * \param values Stores the value of the boundary type + * \param globalPos The global position + */ + void boundaryTypesAtPos(BoundaryTypes &values, + const GlobalPosition &globalPos) const + { + if (globalPos[0] < eps_) + values.setAllDirichlet(); + else + values.setAllNeumann(); + + + // set a dirichlet value for the temperature, use the energy + // equation to set the value + values.setDirichlet(temperatureIdx); + } + + /*! + * \brief Evaluates the boundary conditions for a Dirichlet + * boundary segment + * + * \param values Stores the Dirichlet values for the conservation equations in + * \f$ [ \textnormal{unit of primary variable} ] \f$ + * \param globalPos The global position + */ + void dirichletAtPos(PrimaryVariables &values, const GlobalPosition &globalPos) const + { + // assume a constant density + const Scalar densityW = 1000; + + // assume an intially hydrostatic liquid pressure profile + // note: we subtract rho_w*g*h because g is defined negative + const Scalar pw = 1.0e5 - densityW*this->gravity()[dimWorld-1]*(aquiferDepth_ - globalPos[dimWorld-1]); + + values[Indices::pressureIdx] = pw; + values[saturationIdx] = 0.0; + values[temperatureIdx] = 283.0 + (aquiferDepth_ - globalPos[1])*0.03; + + } + + /*! + * \brief Evaluates the boundary conditions for a Neumann boundary segment. + * + * \param values Stores the Neumann values for the conservation equations in + * \f$ [ \textnormal{unit of conserved quantity} / (m^(dim-1) \cdot s )] \f$ + * \param globalPos The globalPosition of the boundary interface + */ + void neumannAtPos(PrimaryVariables &values, + const GlobalPosition &globalPos) const + { + // initialize values to zero, i.e. no-flow Neumann boundary conditions + values = 0; + values = 0; + + if (this->timeManager().time() + this->timeManager().timeStepSize() < injectionDuration_ + && globalPos[1] < 15 + eps_ && globalPos[1] > 7 - eps_ && globalPos[0] > 0.9*this->bBoxMax()[0]) + { + // inject air. negative values mean injection + values[Indices::contiNEqIdx] = -1e-4; // kg/(s*m^2) + values[Indices::contiWEqIdx] = 0.0; + } + } + + // \} + + + /*! + * \name Volume terms + */ + // \{ + + /*! + * \brief Evaluates the initial values for a control volume + * + * \param values Stores the initial values for the conservation equations in + * \f$ [ \textnormal{unit of primary variables} ] \f$ + * \param globalPos The global position + */ + void initialAtPos(PrimaryVariables &values, const GlobalPosition &globalPos) const + { + // assume a constant density + const Scalar densityW = 1000; + + // assume an intially hydrostatic liquid pressure profile + // note: we subtract rho_w*g*h because g is defined negative + const Scalar pw = 1.0e5 - densityW*this->gravity()[dimWorld-1]*(aquiferDepth_ - globalPos[dimWorld-1]); + + values[Indices::pressureIdx] = pw; + values[saturationIdx] = 0.0; + + values[temperatureIdx] = 283.0 + (aquiferDepth_ - globalPos[1])*0.03; + if (globalPos[0] > 20 - eps_ && globalPos[0] < 30 + eps_ && globalPos[1] > 5 - eps_ && globalPos[1] < 35 + eps_) + values[temperatureIdx] = 380; + } + + +private: + static constexpr Scalar eps_ = 1e-6; + std::string name_; //! Problem name + Scalar aquiferDepth_; //! Depth of the aquifer in m + Scalar injectionDuration_; //! Duration of the injection in seconds +}; +} //end namespace + +#endif diff --git a/tutorial/solution/ex2/injection2p2cproblem.hh b/tutorial/solution/ex2/injection2p2cproblem.hh new file mode 100644 index 0000000000000000000000000000000000000000..716f530614b8e795b2bc030f913feb6e02c5682f --- /dev/null +++ b/tutorial/solution/ex2/injection2p2cproblem.hh @@ -0,0 +1,318 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/***************************************************************************** + * See the file COPYING for full copying permissions. * + * * + * This program is free software: you can redistribute it and/or modify * + * it under the terms of the GNU General Public License as published by * + * the Free Software Foundation, either version 2 of the License, or * + * (at your option) any later version. * + * * + * This program is distributed in the hope that it will be useful, * + * but WITHOUT ANY WARRANTY; without even the implied warranty of * + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * + * GNU General Public License for more details. * + * * + * You should have received a copy of the GNU General Public License * + * along with this program. If not, see . * + *****************************************************************************/ +/*! + * \file + * + * \brief Problem where air is injected under a low permeable layer in a depth of 2700m. + */ +#ifndef DUMUX_INJECTION_2P2C_PROBLEM_HH +#define DUMUX_INJECTION_2P2C_PROBLEM_HH + +#include +#include +#include + +#include "injection2p2cspatialparams.hh" + +// TODO: dumux-course-task +// Include the local residual header +#include "mylocalresidual.hh" + +namespace Dumux +{ + +// foward declaration +template +class Injection2p2cProblem; + +// setup property TypeTag +namespace Properties +{ +// TODO: dumux-course-task +// inherit from MyLocalResidualParams +NEW_TYPE_TAG(Injection2p2cProblem, INHERITS_FROM(TwoPTwoC, InjectionSpatialParams, MyLocalResidualParams)); +NEW_TYPE_TAG(Injection2p2cBoxProblem, INHERITS_FROM(BoxModel, Injection2p2cProblem)); +NEW_TYPE_TAG(Injection2p2pcCCProblem, INHERITS_FROM(CCModel, Injection2p2cProblem)); + +// Set the grid type +SET_TYPE_PROP(Injection2p2cProblem, Grid, Dune::YaspGrid<2>); + +// Set the problem property +SET_TYPE_PROP(Injection2p2cProblem, Problem, Injection2p2cProblem); + +// TODO: dumux-course-task +// change the local residual type to MyTwoPTwoCLocalResidual +SET_TYPE_PROP(Injection2p2cProblem, LocalResidual, MyTwoPTwoCLocalResidual); + +// Set fluid configuration +SET_TYPE_PROP(Injection2p2cProblem, + FluidSystem, + FluidSystems::H2ON2); + +// Define whether mole(true) or mass (false) fractions are used +SET_BOOL_PROP(Injection2p2cProblem, UseMoles, true); + +} // end namespace Properties + +/*! + * \ingroup TwoPTwoCModel + * \ingroup ImplicitTestProblems + * \brief Problem where air is injected under a low permeable layer in a depth of 2700m. + * + * The domain is sized 60m times 40m and consists of two layers, a moderately + * permeable one (\f$ K=10e-12\f$) for \f$ y<22m\f$ and one with a lower permeablility (\f$ K=10e-13\f$) + * in the rest of the domain. + * + * Nitrogen is injected into a water-filled aquifer through a well. First, we inject for one month. + * Then, we continue simulating the development of the nitrogen plume for 10 years. + * This is realized with a Neumann boundary condition at the right boundary + * (\f$ 7m +class Injection2p2cProblem : public ImplicitPorousMediaProblem +{ + typedef ImplicitPorousMediaProblem ParentType; + + typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView; + typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar; + typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem; + + // grid world dimension + static constexpr auto dimWorld = GridView::dimensionworld; + + typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices; + typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables; + typedef typename GET_PROP_TYPE(TypeTag, BoundaryTypes) BoundaryTypes; + typedef typename GET_PROP_TYPE(TypeTag, TimeManager) TimeManager; + typedef typename GET_PROP_TYPE(TypeTag, MaterialLaw) MaterialLaw; + typedef Dune::FieldVector GlobalPosition; + +public: + /*! + * \brief The constructor + * + * \param timeManager The time manager + * \param gridView The grid view + */ + Injection2p2cProblem(TimeManager &timeManager, const GridView &gridView) + : ParentType(timeManager, gridView) + { + // initialize the tables of the fluid system + FluidSystem::init(/*tempMin=*/273.15, + /*tempMax=*/423.15, + /*numTemp=*/50, + /*pMin=*/0.0, + /*pMax=*/30e6, + /*numP=*/300); + + // name of the problem and output file + name_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, std::string, Problem, Name); + // depth of the aquifer, units: m + aquiferDepth_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Problem, AquiferDepth); + // inflow rate of nitrogen water vapor mixture, units: kg/(s m^2) + totalAreaSpecificInflow_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Problem, TotalAreaSpecificInflow); + // the duration of the injection, units: second + injectionDuration_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Problem, InjectionDuration); + } + + /*! + * \brief User defined output after the time integration + * + * Will be called diretly after the time integration + */ + void postTimeStep() + { + // Calculate storage terms + PrimaryVariables storageW, storageN; + this->model().globalPhaseStorage(storageW, Indices::wPhaseIdx); + this->model().globalPhaseStorage(storageN, Indices::nPhaseIdx); + + // Write mass balance information for rank 0 + if (this->gridView().comm().rank() == 0) + { + std::cout <<"Storage: wetting=[" << storageW << "]" + << " nonwetting=[" << storageN << "]" << std::endl; + } + } + + /*! + * \name Problem parameters + */ + // \{ + + /*! + * \brief Returns the problem name + * + * This is used as a prefix for files generated by the simulation. + */ + const std::string& name() const + { return name_; } + + /*! + * \brief Returns the temperature in \f$ K \f$ + */ + Scalar temperature() const + { return 273.15 + 30; } + + /*! + * \brief Returns the source term + * + * \param values Stores the source values for the conservation equations in + * \f$ [ \textnormal{unit of primary variable} / (m^\textrm{dim} \cdot s )] \f$ + * \param globalPos The global position + */ + void sourceAtPos(PrimaryVariables &values, + const GlobalPosition &globalPos) const + { values = 0; } + + // \} + + /*! + * \name Boundary conditions + */ + // \{ + + /*! + * \brief Specifies which kind of boundary condition should be + * used for which equation on a given boundary segment + * + * \param values Stores the value of the boundary type + * \param globalPos The global position + */ + void boundaryTypesAtPos(BoundaryTypes &values, + const GlobalPosition &globalPos) const + { + // Set Dirichlet at the bottom of the domain + if (globalPos[dimWorld-1] < eps_) + { + values.setAllDirichlet(); + } + + // and Neuman boundary conditions everywhere else + // note that we don't differentiate between Neumann and Robin boundary types + else + { + values.setAllNeumann(); + } + } + + /*! + * \brief Evaluates the boundary conditions for a Dirichlet + * boundary segment + * + * \param values Stores the Dirichlet values for the conservation equations in + * \f$ [ \textnormal{unit of primary variable} ] \f$ + * \param globalPos The global position + */ + void dirichletAtPos(PrimaryVariables &values, const GlobalPosition &globalPos) const + { initialAtPos(values, globalPos); } + + /*! + * \brief Evaluates the boundary conditions for a Neumann boundary segment. + * + * \param values Stores the Neumann values for the conservation equations in + * \f$ [ \textnormal{unit of conserved quantity} / (m^(dim-1) \cdot s )] \f$ + * \param globalPos The globalPosition of the boundary interface + */ + void neumannAtPos(PrimaryVariables &values, + const GlobalPosition &globalPos) const + { + // initialize values to zero, i.e. no-flow Neumann boundary conditions + values = 0; + + //if we are inside the injection zone set inflow Neumann boundary conditions + if (this->timeManager().time() + this->timeManager().timeStepSize() < injectionDuration_ + && globalPos[1] < 15 + eps_ && globalPos[1] > 7 - eps_ && globalPos[0] > 0.9*this->bBoxMax()[0]) + { + // set the Neumann values for the Nitrogen component balance + // convert from units kg/(s*m^2) to mole/(s*m^2) + values[Indices::contiNEqIdx] = -totalAreaSpecificInflow_/FluidSystem::molarMass(FluidSystem::nCompIdx); + values[Indices::contiWEqIdx] = 0.0; + } + } + + // \} + + /*! + * \name Volume terms + */ + // \{ + + /*! + * \brief Evaluates the initial values for a control volume + * + * \param values Stores the initial values for the conservation equations in + * \f$ [ \textnormal{unit of primary variables} ] \f$ + * \param globalPos The global position + */ + void initialAtPos(PrimaryVariables &values, const GlobalPosition &globalPos) const + { + // get the water density at atmospheric conditions + const Scalar densityW = FluidSystem::H2O::liquidDensity(temperature(), 1.0e5); + + // assume an intially hydrostatic liquid pressure profile + // note: we subtract rho_w*g*h because g is defined negative + const Scalar pw = 1.0e5 - densityW*this->gravity()[dimWorld-1]*(aquiferDepth_ - globalPos[dimWorld-1]); + + // initially we have some nitrogen dissolved + // saturation mole fraction would be + // moleFracLiquidN2 = (pw + pc + p_vap^sat)/henry; + const Scalar moleFracLiquidN2 = pw*0.95/BinaryCoeff::H2O_N2::henry(temperature()); + + // note that because we start with a single phase system the primary variables + // are pl and x^w_N2. This will switch as soon after we start injecting to a two + // phase system so the primary variables will be pl and Sn (non-wetting saturation). + values[Indices::switchIdx] = moleFracLiquidN2; + values[Indices::pressureIdx] = pw; + } + + /*! + * \brief Return the initial phase state inside a control volume. + * + * \param globalPos The global position + * \note we start with a single phase system + */ + int initialPhasePresenceAtPos(const GlobalPosition &globalPos) const + { return Indices::wPhaseOnly; } + + // \} + + //! If we should write restart files + bool shouldWriteRestartFile() const + { return false; } + +private: + static constexpr Scalar eps_ = 1e-6; + std::string name_; //! Problem name + Scalar aquiferDepth_; //! Depth of the aquifer in m + Scalar totalAreaSpecificInflow_; //! Area specific inflow rate in mole/(s*m^2) + Scalar injectionDuration_; //! Duration of the injection in seconds +}; + +} // end namespace Dumux + +#endif diff --git a/tutorial/solution/ex2/injection2p2cspatialparams.hh b/tutorial/solution/ex2/injection2p2cspatialparams.hh new file mode 100644 index 0000000000000000000000000000000000000000..de2eb158567b8dcc00500b18aba9f96adf096764 --- /dev/null +++ b/tutorial/solution/ex2/injection2p2cspatialparams.hh @@ -0,0 +1,223 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/***************************************************************************** + * See the file COPYING for full copying permissions. * + * * + * This program is free software: you can redistribute it and/or modify * + * it under the terms of the GNU General Public License as published by * + * the Free Software Foundation, either version 2 of the License, or * + * (at your option) any later version. * + * * + * This program is distributed in the hope that it will be useful, * + * but WITHOUT ANY WARRANTY; without even the implied warranty of * + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * + * GNU General Public License for more details. * + * * + * You should have received a copy of the GNU General Public License * + * along with this program. If not, see . * + *****************************************************************************/ +/*! + * \file + * + * \brief Definition of the spatial parameters for the injection problem + * which uses the isothermal two-phase two-component + * fully implicit model. + */ + +#ifndef DUMUX_INJECTION_SPATIAL_PARAMS_HH +#define DUMUX_INJECTION_SPATIAL_PARAMS_HH + +#include +#include +#include +// TODO: dumux-course-task +// Inlcude your own material law +#include "mymateriallaw.hh" + +#include +#include + +#include + +namespace Dumux +{ + +// forward declaration +template +class InjectionSpatialParams; + +// setup property TypeTag +namespace Properties +{ +// The spatial parameters TypeTag +NEW_TYPE_TAG(InjectionSpatialParams); + +// Set the spatial parameters +SET_TYPE_PROP(InjectionSpatialParams, SpatialParams, InjectionSpatialParams); + +// TODO: dumux-course-task +// Use your own material law instead +// Set the material law parameterized by absolute saturations +SET_PROP(InjectionSpatialParams, MaterialLaw) +{ + using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar); + // using type = EffToAbsLaw>; + using type = EffToAbsLaw>; +}; + +} // end namespace Properties + +/*! + * \ingroup TwoPTwoCModel + * \ingroup ImplicitTestProblems + * \brief Definition of the spatial parameters for the injection problem + * which uses the isothermal two-phase two-component + * fully implicit model. + */ +template +class InjectionSpatialParams : public ImplicitSpatialParams +{ + typedef ImplicitSpatialParams ParentType; + typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView; + typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar; + typedef typename GridView::ctype CoordinateType; + + enum { + dim=GridView::dimension, + dimWorld=GridView::dimensionworld + }; + + typedef Dune::FieldVector GlobalPosition; + typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry; + typedef typename GridView::template Codim<0>::Entity Element; + typedef typename GET_PROP_TYPE(TypeTag, MaterialLaw) MaterialLaw; + typedef typename MaterialLaw::Params MaterialLawParams; + +public: + + /*! + * \brief The constructor + * + * \param gridView The grid view + */ + InjectionSpatialParams(const GridView &gridView) + : ParentType(gridView) + { + aquiferHeightFromBottom_ = 30.0; + + // intrinsic permeabilities + aquitardK_ = GET_RUNTIME_PARAM(TypeTag, Scalar, SpatialParams.PermeabilityAquitard); + aquiferK_ = GET_RUNTIME_PARAM(TypeTag, Scalar, SpatialParams.PermeabilityAquifer); + + // porosities + aquitardPorosity_ = 0.2; + aquiferPorosity_ = 0.4; + + // residual saturations + aquitardMaterialParams_.setSwr(0.2); + aquitardMaterialParams_.setSnr(0.0); + aquiferMaterialParams_.setSwr(0.2); + aquiferMaterialParams_.setSnr(0.0); + + // parameters for the Brooks-Corey law + aquitardMaterialParams_.setPe(GET_RUNTIME_PARAM(TypeTag, Scalar, SpatialParams.EntryPressureAquitard)); + aquiferMaterialParams_.setPe(GET_RUNTIME_PARAM(TypeTag, Scalar, SpatialParams.EntryPressureAquifer)); + aquitardMaterialParams_.setLambda(2.0); + aquiferMaterialParams_.setLambda(2.0); + + // plot the material laws using gnuplot and exit + if (GET_RUNTIME_PARAM(TypeTag, bool, Problem.OnlyPlotMaterialLaws)) + { + plotMaterialLaws(); + exit(0); + } + } + + /*! + * \brief Returns the intrinsic permeability tensor \f$[m^2]\f$ + * + * \param element The finite element + * \param fvGeometry The finite volume geometry of the element + * \param scvIdx The local index of the sub-control volume + */ + Scalar intrinsicPermeability(const Element &element, + const FVElementGeometry &fvGeometry, + const int scvIdx) const + { + const GlobalPosition &globalPos = fvGeometry.subContVol[scvIdx].global; + if (isInAquitard_(globalPos)) + return aquitardK_; + return aquiferK_; + } + + /*! + * \brief Returns the porosity \f$[-]\f$ + * + * \param element The finite element + * \param fvGeometry The finite volume geometry of the element + * \param scvIdx The local index of the sub-control volume + */ + Scalar porosity(const Element &element, + const FVElementGeometry &fvGeometry, + const int scvIdx) const + { + const GlobalPosition &globalPos = fvGeometry.subContVol[scvIdx].global; + if (isInAquitard_(globalPos)) + return aquitardPorosity_; + return aquiferPorosity_; + } + + + /*! + * \brief Returns the parameter object for the capillary-pressure/ + * saturation material law + * + * \param element The finite element + * \param fvGeometry The finite volume geometry of the element + * \param scvIdx The local index of the sub-control volume + */ + const MaterialLawParams& materialLawParams(const Element &element, + const FVElementGeometry &fvGeometry, + const int scvIdx) const + { + const GlobalPosition &globalPos = fvGeometry.subContVol[scvIdx].global; + if (isInAquitard_(globalPos)) + return aquitardMaterialParams_; + return aquiferMaterialParams_; + } + + /*! + * \brief Creates a gnuplot output of the pc-Sw curve + */ + void plotMaterialLaws() + { + PlotMaterialLaw plotMaterialLaw; + GnuplotInterface gnuplot; + plotMaterialLaw.addpcswcurve(gnuplot, aquitardMaterialParams_, 0.2, 1.0, "upper layer (fine, aquitard)", "w lp"); + plotMaterialLaw.addpcswcurve(gnuplot, aquiferMaterialParams_, 0.2, 1.0, "lower layer (coarse, aquifer)", "w l"); + gnuplot.setOption("set xrange [0:1]"); + gnuplot.setOption("set label \"residual\\nsaturation\" at 0.1,100000 center"); + gnuplot.plot("pc-Sw"); + } + +private: + + static constexpr Scalar eps_ = 1e-6; + + bool isInAquitard_(const GlobalPosition &globalPos) const + { return globalPos[dimWorld-1] > aquiferHeightFromBottom_ + eps_; } + + Scalar aquitardK_; + Scalar aquiferK_; + Scalar aquiferHeightFromBottom_; + + Scalar aquitardPorosity_; + Scalar aquiferPorosity_; + + MaterialLawParams aquitardMaterialParams_; + MaterialLawParams aquiferMaterialParams_; +}; + +} // end namespace Dumux + +#endif diff --git a/tutorial/solution/ex2/mylocalresidual.hh b/tutorial/solution/ex2/mylocalresidual.hh new file mode 100644 index 0000000000000000000000000000000000000000..8995eb80d5c6c14898caabdc6afd99af60011e01 --- /dev/null +++ b/tutorial/solution/ex2/mylocalresidual.hh @@ -0,0 +1,545 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/***************************************************************************** + * See the file COPYING for full copying permissions. * + * * + * This program is free software: you can redistribute it and/or modify * + * it under the terms of the GNU General Public License as published by * + * the Free Software Foundation, either version 2 of the License, or * + * (at your option) any later version. * + * * + * This program is distributed in the hope that it will be useful, * + * but WITHOUT ANY WARRANTY; without even the implied warranty of * + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * + * GNU General Public License for more details. * + * * + * You should have received a copy of the GNU General Public License * + * along with this program. If not, see . * + *****************************************************************************/ +/*! + * \file + * + * \brief Element-wise calculation of the Jacobian matrix for problems + * using the two-phase two-component fully implicit model. + */ + +#ifndef DUMUX_MY_2P2C_LOCAL_RESIDUAL_HH +#define DUMUX_MY_2P2C_LOCAL_RESIDUAL_HH + +#include + +namespace Dumux +{ + +// TODO: dumux-course-task +// implement new TypeTag MyLocalResidualParams +namespace Properties +{ +NEW_TYPE_TAG(MyLocalResidualParams); +NEW_PROP_TAG(ProblemEnableDiffusion); +SET_BOOL_PROP(MyLocalResidualParams, ProblemEnableDiffusion, true); +} // end namespace Properties + +/*! + * \ingroup TwoPTwoCModel + * \ingroup ImplicitLocalResidual + * \brief Element-wise calculation of the Jacobian matrix for problems + * using the two-phase two-component fully implicit model. + * + * This class is used to fill the gaps in ImplicitLocalResidual for the + * two-phase two-component flow. + */ +template +class MyTwoPTwoCLocalResidual: public GET_PROP_TYPE(TypeTag, BaseLocalResidual) +{ + protected: + typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar; + typedef typename GET_PROP_TYPE(TypeTag, LocalResidual) Implementation; + typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem; + typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables; + typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry; + typedef typename GET_PROP_TYPE(TypeTag, VolumeVariables) VolumeVariables; + typedef typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables) ElementVolumeVariables; + typedef typename GET_PROP_TYPE(TypeTag, ElementBoundaryTypes) ElementBoundaryTypes; + typedef typename GET_PROP_TYPE(TypeTag, FluxVariables) FluxVariables; + enum + { + numPhases = GET_PROP_VALUE(TypeTag, NumPhases), + numComponents = GET_PROP_VALUE(TypeTag, NumComponents) + }; + + typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices; + enum + { + contiWEqIdx = Indices::contiWEqIdx, + contiNEqIdx = Indices::contiNEqIdx, + wPhaseIdx = Indices::wPhaseIdx, + nPhaseIdx = Indices::nPhaseIdx, + wCompIdx = Indices::wCompIdx, + nCompIdx = Indices::nCompIdx + }; + + typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView; + typedef typename GridView::template Codim<0>::Entity Element; + + static constexpr unsigned int replaceCompEqIdx = GET_PROP_VALUE(TypeTag, ReplaceCompEqIdx); + + //! Property that defines whether mole or mass fractions are used + static const bool useMoles = GET_PROP_VALUE(TypeTag, UseMoles); + + public: + /*! + * \brief Constructor + * + * Sets the mass upwind weight. + */ + MyTwoPTwoCLocalResidual() + { + // retrieve the upwind weight for the mass conservation equations. Use the value + // specified via the property system as default, and overwrite + // it by the run-time parameter from the Dune::ParameterTree + massUpwindWeight_ = GET_PARAM_FROM_GROUP(TypeTag, Scalar, Implicit, MassUpwindWeight); + + // TODO: dumux-course-task + // get parameter Problem.EnableDiffusion + enableDiffusion_ = GET_PARAM_FROM_GROUP(TypeTag, bool, Problem, EnableDiffusion); + } + + /*! + * \brief Evaluate the storage term of the current solution in a + * single phase. + * + * \param element The element + * \param phaseIdx The index of the fluid phase + */ + void evalPhaseStorage(const Element &element, const int phaseIdx) + { + FVElementGeometry fvGeometry; + fvGeometry.update(this->gridView_(), element); + ElementBoundaryTypes bcTypes; + bcTypes.update(this->problem_(), element, fvGeometry); + ElementVolumeVariables elemVolVars; + elemVolVars.update(this->problem_(), element, fvGeometry, false); + + this->storageTerm_.resize(fvGeometry.numScv); + this->storageTerm_ = 0; + + this->elemPtr_ = &element; + this->fvElemGeomPtr_ = &fvGeometry; + this->bcTypesPtr_ = &bcTypes; + this->prevVolVarsPtr_ = 0; + this->curVolVarsPtr_ = &elemVolVars; + evalPhaseStorage_(phaseIdx); + } + + /*! + * \brief Evaluate the amount of all conservation quantities + * (e.g. phase mass) within a sub-control volume. + * + * \param storage The mass of the component within the sub-control volume + * \param scvIdx The sub-control-volume index + * \param usePrevSol Based on usePrevSol solution of current or previous time step is used + * + * The result should be averaged over the volume (e.g. phase mass + * inside a sub-control volume divided by the volume) + */ + void computeStorage(PrimaryVariables &storage, const int scvIdx, bool usePrevSol) const + { + // if flag usePrevSol is set, the solution from the previous + // time step is used, otherwise the current solution is + // used. The secondary variables are used accordingly. This + // is required to compute the derivative of the storage term + // using the implicit Euler method. + const ElementVolumeVariables &elemVolVars = usePrevSol ? this->prevVolVars_() + : this->curVolVars_(); + const VolumeVariables &volVars = elemVolVars[scvIdx]; + + // compute storage term of all components within all phases + storage = 0; + if(useMoles) // mole-fraction formulation + { + for (unsigned int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) + { + for (unsigned int compIdx = contiCompIdx1_(); compIdx <= contiCompIdx2_(); ++compIdx) + { + unsigned int eqIdx = (compIdx == wCompIdx) ? contiWEqIdx : contiNEqIdx; + storage[eqIdx] += volVars.molarDensity(phaseIdx) + * volVars.saturation(phaseIdx) + * volVars.moleFraction(phaseIdx, compIdx); + } + // this is only processed if one component mass balance equation + // is replaced by the total mass balance equation + if (replaceCompEqIdx < numComponents) + storage[replaceCompEqIdx] += + volVars.molarDensity(phaseIdx) + * volVars.saturation(phaseIdx); + } + storage *= volVars.porosity(); + } + else // mass-fraction formulation + { + for (unsigned int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) + { + for (unsigned int compIdx = contiCompIdx1_(); compIdx <= contiCompIdx2_(); ++compIdx) + { + unsigned int eqIdx = (compIdx == wCompIdx) ? contiWEqIdx : contiNEqIdx; + storage[eqIdx] += volVars.density(phaseIdx) + * volVars.saturation(phaseIdx) + * volVars.massFraction(phaseIdx, compIdx); + } + // this is only processed if one component mass balance equation + // is replaced by the total mass balance equation + if (replaceCompEqIdx < numComponents) + storage[replaceCompEqIdx] += + volVars.density(phaseIdx) + * volVars.saturation(phaseIdx); + } + storage *= volVars.porosity(); + } + } + + /*! + * \brief Evaluates the total flux of all conservation quantities + * over a face of a sub-control volume. + * + * \param flux The flux over the sub-control-volume face for each component + * \param fIdx The index of the sub-control-volume face + * \param onBoundary Evaluate flux at inner sub-control-volume face or on a boundary face + */ + void computeFlux(PrimaryVariables &flux, const int fIdx, bool onBoundary=false) const + { + // update the flux variables + FluxVariables fluxVars; + fluxVars.update(this->problem_(), + this->element_(), + this->fvGeometry_(), + fIdx, + this->curVolVars_(), + onBoundary); + + flux = 0; + + // add advective fluxes + computeAdvectiveFlux(flux, fluxVars); + + // add diffusive fluxes + // TODO: dumux-course-task + // only add diffusive fluxes if diffusion is enabled + if (enableDiffusion_) + computeDiffusiveFlux(flux, fluxVars); + } + + /*! + * \brief Evaluates the advective mass flux of all components over + * a face of a sub-control volume. + * + * \param flux The flux over the sub-control-volume face for each component + * \param fluxVars The flux variables at the current sub-control-volume face + */ + void computeAdvectiveFlux(PrimaryVariables &flux, const FluxVariables &fluxVars) const + { + //////// + // advective fluxes of all components in all phases + //////// + + if(useMoles) // mole-fraction formulation + { + for (unsigned int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) + { + // data attached to upstream and the downstream vertices + // of the current phase + const VolumeVariables &up = + this->curVolVars_(fluxVars.upstreamIdx(phaseIdx)); + const VolumeVariables &dn = + this->curVolVars_(fluxVars.downstreamIdx(phaseIdx)); + + for (unsigned int compIdx = contiCompIdx1_(); compIdx <= contiCompIdx2_(); ++compIdx) + { + unsigned int eqIdx = (compIdx == wCompIdx) ? contiWEqIdx : contiNEqIdx; + // add advective flux of current component in current + // phase + if (massUpwindWeight_ > 0.0) + // upstream vertex + flux[eqIdx] += + fluxVars.volumeFlux(phaseIdx) + * massUpwindWeight_ + * up.molarDensity(phaseIdx) + * up.moleFraction(phaseIdx, compIdx); + if (massUpwindWeight_ < 1.0) + // downstream vertex + flux[eqIdx] += + fluxVars.volumeFlux(phaseIdx) + * (1 - massUpwindWeight_) + * dn.molarDensity(phaseIdx) + * dn.moleFraction(phaseIdx, compIdx); + + Valgrind::CheckDefined(fluxVars.volumeFlux(phaseIdx)); + Valgrind::CheckDefined(up.molarDensity(phaseIdx)); + Valgrind::CheckDefined(up.moleFraction(phaseIdx, compIdx)); + Valgrind::CheckDefined(dn.molarDensity(phaseIdx)); + Valgrind::CheckDefined(dn.moleFraction(phaseIdx, compIdx)); + } + // flux of the total mass balance; + // this is only processed if one component mass balance equation + // is replaced by a total mass balance equation + if (replaceCompEqIdx < numComponents) + { + // upstream vertex + if (massUpwindWeight_ > 0.0) + flux[replaceCompEqIdx] += + fluxVars.volumeFlux(phaseIdx) + * massUpwindWeight_ + * up.molarDensity(phaseIdx); + // downstream vertex + if (massUpwindWeight_ < 1.0) + flux[replaceCompEqIdx] += + fluxVars.volumeFlux(phaseIdx) + * (1 - massUpwindWeight_) + * dn.molarDensity(phaseIdx); + Valgrind::CheckDefined(fluxVars.volumeFlux(phaseIdx)); + Valgrind::CheckDefined(up.molarDensity(phaseIdx)); + Valgrind::CheckDefined(dn.molarDensity(phaseIdx)); + + } + + } + } + else // mass-fraction formulation + { + for (unsigned int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) + { + // data attached to upstream and downstream vertices + // of the current phase + const VolumeVariables &up = + this->curVolVars_(fluxVars.upstreamIdx(phaseIdx)); + const VolumeVariables &dn = + this->curVolVars_(fluxVars.downstreamIdx(phaseIdx)); + + for (unsigned int compIdx = contiCompIdx1_(); compIdx <= contiCompIdx2_(); ++compIdx) + { + unsigned int eqIdx = (compIdx == wCompIdx) ? contiWEqIdx : contiNEqIdx; + // add advective flux of current component in current + // phase + if (massUpwindWeight_ > 0.0) + // upstream vertex + flux[eqIdx] += + fluxVars.volumeFlux(phaseIdx) + * massUpwindWeight_ + * up.density(phaseIdx) + * up.massFraction(phaseIdx, compIdx); + if (massUpwindWeight_ < 1.0) + // downstream vertex + flux[eqIdx] += + fluxVars.volumeFlux(phaseIdx) + * (1 - massUpwindWeight_) + * dn.density(phaseIdx) + * dn.massFraction(phaseIdx, compIdx); + + Valgrind::CheckDefined(fluxVars.volumeFlux(phaseIdx)); + Valgrind::CheckDefined(up.density(phaseIdx)); + Valgrind::CheckDefined(up.massFraction(phaseIdx, compIdx)); + Valgrind::CheckDefined(dn.density(phaseIdx)); + Valgrind::CheckDefined(dn.massFraction(phaseIdx, compIdx)); + } + // flux of the total mass balance; + // this is only processed if one component mass balance equation + // is replaced by a total mass balance equation + if (replaceCompEqIdx < numComponents) + { + // upstream vertex + if (massUpwindWeight_ > 0.0) + flux[replaceCompEqIdx] += + fluxVars.volumeFlux(phaseIdx) + * massUpwindWeight_ + * up.density(phaseIdx); + // downstream vertex + if (massUpwindWeight_ < 1.0) + flux[replaceCompEqIdx] += + fluxVars.volumeFlux(phaseIdx) + * (1 - massUpwindWeight_) + * dn.density(phaseIdx); + Valgrind::CheckDefined(fluxVars.volumeFlux(phaseIdx)); + Valgrind::CheckDefined(up.density(phaseIdx)); + Valgrind::CheckDefined(dn.density(phaseIdx)); + + } + + } + } + } + + /*! + * \brief Evaluates the diffusive mass flux of all components over + * a face of a sub-control volume. + * + * \param flux The flux over the sub-control-volume face for each component + * \param fluxVars The flux variables at the current sub-control-volume face + */ + void computeDiffusiveFlux(PrimaryVariables &flux, const FluxVariables &fluxVars) const + + { + if(useMoles) // mole-fraction formulation + { + // add diffusive flux of gas component in liquid phase + Scalar tmp = -(fluxVars.moleFractionGrad(wPhaseIdx)*fluxVars.face().normal) + * fluxVars.porousDiffCoeff(wPhaseIdx) + * fluxVars.molarDensity(wPhaseIdx); + // add the diffusive fluxes only to the component mass balance + if (replaceCompEqIdx != contiNEqIdx) + flux[contiNEqIdx] += tmp; + if (replaceCompEqIdx != contiWEqIdx) + flux[contiWEqIdx] -= tmp; + + // add diffusive flux of liquid component in non-wetting phase + tmp = -(fluxVars.moleFractionGrad(nPhaseIdx)*fluxVars.face().normal) + * fluxVars.porousDiffCoeff(nPhaseIdx) + * fluxVars.molarDensity(nPhaseIdx); + // add the diffusive fluxes only to the component mass balance + if (replaceCompEqIdx != contiWEqIdx) + flux[contiWEqIdx] += tmp; + if (replaceCompEqIdx != contiNEqIdx) + flux[contiNEqIdx] -= tmp; + } + else // mass-fraction formulation + { + // add diffusive flux of gas component in wetting phase + Scalar tmp = -(fluxVars.moleFractionGrad(wPhaseIdx)*fluxVars.face().normal) + * fluxVars.porousDiffCoeff(wPhaseIdx) + * fluxVars.molarDensity(wPhaseIdx); + // add the diffusive fluxes to the component mass balance and total mass balance + if (replaceCompEqIdx < numComponents) + { + flux[replaceCompEqIdx] += tmp * FluidSystem::molarMass(nCompIdx); + flux[replaceCompEqIdx] -= tmp * FluidSystem::molarMass(wCompIdx); + } + if (replaceCompEqIdx != contiWEqIdx) + { + flux[contiWEqIdx] -= tmp * FluidSystem::molarMass(wCompIdx); + } + if (replaceCompEqIdx != contiNEqIdx) + { + flux[contiNEqIdx] += tmp * FluidSystem::molarMass(nCompIdx); + } + + // add diffusive fluxes of liquid component in non-wetting phase + tmp = -(fluxVars.moleFractionGrad(nPhaseIdx)*fluxVars.face().normal) + * fluxVars.porousDiffCoeff(nPhaseIdx) + * fluxVars.molarDensity(nPhaseIdx); + // add the diffusive fluxes to the component mass balance and total mass balance + if (replaceCompEqIdx < numComponents) + { + flux[replaceCompEqIdx] += tmp * FluidSystem::molarMass(wCompIdx); + flux[replaceCompEqIdx] -= tmp * FluidSystem::molarMass(nCompIdx); + } + if (replaceCompEqIdx != contiWEqIdx) + { + flux[contiWEqIdx] += tmp * FluidSystem::molarMass(wCompIdx); + } + if (replaceCompEqIdx != contiNEqIdx) + { + flux[contiNEqIdx] -= tmp * FluidSystem::molarMass(nCompIdx); + } + } + } + + protected: + void evalPhaseStorage_(const int phaseIdx) + { + if(useMoles) // mole-fraction formulation + { + // evaluate the storage terms of a single phase + for (int i=0; i < this->fvGeometry_().numScv; i++) { + PrimaryVariables &storage = this->storageTerm_[i]; + const ElementVolumeVariables &elemVolVars = this->curVolVars_(); + const VolumeVariables &volVars = elemVolVars[i]; + + // compute storage term of all components within all phases + storage = 0; + for (int compIdx = 0; compIdx < numComponents; ++compIdx) + { + int eqIdx = (compIdx == wCompIdx) ? contiWEqIdx : contiNEqIdx; + storage[eqIdx] += volVars.molarDensity(phaseIdx) + * volVars.saturation(phaseIdx) + * volVars.moleFraction(phaseIdx, compIdx); + } + + storage *= volVars.porosity(); + storage *= this->fvGeometry_().subContVol[i].volume; + } + } + else // mass-fraction formulation + { + // evaluate the storage terms of a single phase + for (int i=0; i < this->fvGeometry_().numScv; i++) { + PrimaryVariables &storage = this->storageTerm_[i]; + const ElementVolumeVariables &elemVolVars = this->curVolVars_(); + const VolumeVariables &volVars = elemVolVars[i]; + + // compute storage term of all components within all phases + storage = 0; + for (int compIdx = 0; compIdx < numComponents; ++compIdx) + { + int eqIdx = (compIdx == wCompIdx) ? contiWEqIdx : contiNEqIdx; + storage[eqIdx] += volVars.density(phaseIdx) + * volVars.saturation(phaseIdx) + * volVars.massFraction(phaseIdx, compIdx); + } + + storage *= volVars.porosity(); + storage *= this->fvGeometry_().subContVol[i].volume; + } + } + } + + /*! + * \brief Returns the equation index of the first mass-balance equation + * of the component (used for loops) + * + * Returns the equation index of the first mass-balance equation + * of the component (used for loops) if one component mass balance + * is replaced by the total mass balance, this is the index + * of the remaining component mass-balance equation. + */ + unsigned int contiCompIdx1_() const { + switch (replaceCompEqIdx) + { + case contiWEqIdx: return contiNEqIdx; + case contiNEqIdx: return contiWEqIdx; + default: return 0; + } + } + + /*! + * \brief Returns the equation index of the second mass balance + * of the component (used for loops) + * + * Returns the equation index of the second mass balance + * of the component (used for loops) + * if one component mass balance is replaced by the total mass balance + * (replaceCompEqIdx < 2), this index is the same as contiCompIdx1(). + */ + unsigned int contiCompIdx2_() const { + switch (replaceCompEqIdx) + { + case contiWEqIdx: return contiNEqIdx; + case contiNEqIdx: return contiWEqIdx; + default: return numComponents-1; + } + } + + Implementation *asImp_() + { return static_cast (this); } + const Implementation *asImp_() const + { return static_cast (this); } + + private: + Scalar massUpwindWeight_; + + // TODO: dumux-course-task + // add private member enableDiffusion_ + bool enableDiffusion_; +}; + +} // end namespace Dumux + +#endif diff --git a/tutorial/solution/ex2/mymateriallaw.hh b/tutorial/solution/ex2/mymateriallaw.hh new file mode 100644 index 0000000000000000000000000000000000000000..dbd58c8ff9b29fc2f952112548eaf5ed78a297dd --- /dev/null +++ b/tutorial/solution/ex2/mymateriallaw.hh @@ -0,0 +1,115 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/***************************************************************************** + * See the file COPYING for full copying permissions. * + * * + * This program is free software: you can redistribute it and/or modify * + * it under the terms of the GNU General Public License as published by * + * the Free Software Foundation, either version 2 of the License, or * + * (at your option) any later version. * + * * + * This program is distributed in the hope that it will be useful, * + * but WITHOUT ANY WARRANTY; without even the implied warranty of * + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * + * GNU General Public License for more details. * + * * + * You should have received a copy of the GNU General Public License * + * along with this program. If not, see . * + *****************************************************************************/ +/*! + * \file + * + * \brief Implementation of the capillary pressure and + * relative permeability <-> saturation relations. + * + */ +#ifndef DUMUX_MY_MATERIAL_LAW_HH +#define DUMUX_MY_MATERIAL_LAW_HH + +#include +#include + +namespace Dumux +{ +/*! + * \ingroup fluidmatrixinteractionslaws + * \note a simple material law using the BrooksCoreyParams + */ +template > +class MyMaterialLaw +{ +public: + typedef ParamsT Params; + typedef typename Params::Scalar Scalar; + + /*! + * \brief The capillary pressure-saturation curve + * \param swe saturation of the wetting phase \f$\mathrm{[\overline{S}_w]}\f$ + * \param params A container object that is populated with the appropriate coefficients for the respective law. + * Therefore, in the (problem specific) spatialParameters first, the material law is chosen, + and then the params container is constructed accordingly. Afterwards the values are set there, too. + * \return capillary pressure + * TODO: dumux-course-task + * Implement the pc(swe) function + */ + static Scalar pc(const Params ¶ms, Scalar swe) + { + return 1.0e5*(1.0-swe) + params.pe(); + } + + /*! + * \brief The relative permeability for the wetting phase of + * the medium implied by the Brooks-Corey + * parameterization. + * + * \param swe The mobile saturation of the wetting phase. + * \param params A container object that is populated with the appropriate coefficients for the respective law. + * Therefore, in the (problem specific) spatialParameters first, the material law is chosen, + * and then the params container is constructed accordingly. Afterwards the values are set there, too. + * \return Relative permeability of the wetting phase calculated as implied by Brooks & Corey. + * + * \note Instead of undefined behaviour if pc is not in the valid range, we return a valid number, + * by clamping the input. + */ + static Scalar krw(const Params ¶ms, Scalar swe) + { + using std::pow; + using std::min; + using std::max; + + swe = min(max(swe, 0.0), 1.0); // the equation below is only defined for 0.0 <= swe <= 1.0 + + return pow(swe, 2.0/params.lambda() + 3); + } + + /*! + * \brief The relative permeability for the non-wetting phase of + * the medium as implied by the Brooks-Corey + * parameterization. + * + * \param swe The mobile saturation of the wetting phase. + * \param params A container object that is populated with the appropriate coefficients for the respective law. + * Therefore, in the (problem specific) spatialParameters first, the material law is chosen, and then the params container + * is constructed accordingly. Afterwards the values are set there, too. + * \return Relative permeability of the non-wetting phase calculated as implied by Brooks & Corey. + * + * \note Instead of undefined behaviour if pc is not in the valid range, we return a valid number, + * by clamping the input. + */ + static Scalar krn(const Params ¶ms, Scalar swe) + { + using std::pow; + using std::min; + using std::max; + + swe = min(max(swe, 0.0), 1.0); // the equation below is only defined for 0.0 <= swe <= 1.0 + + const Scalar exponent = 2.0/params.lambda() + 1; + const Scalar tmp = 1.0 - swe; + return tmp*tmp*(1.0 - pow(swe, exponent)); + } +}; + +} // end namespace Dumux + +#endif diff --git a/tutorial/solution/ex3/h2omycompressiblecomponent.hh b/tutorial/solution/ex3/h2omycompressiblecomponent.hh new file mode 100644 index 0000000000000000000000000000000000000000..47822a8647b1d5280379229de1456ff43258b99d --- /dev/null +++ b/tutorial/solution/ex3/h2omycompressiblecomponent.hh @@ -0,0 +1,455 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/***************************************************************************** + * See the file COPYING for full copying permissions. * + * * + * This program is free software: you can redistribute it and/or modify * + * it under the terms of the GNU General Public License as published by * + * the Free Software Foundation, either version 2 of the License, or * + * (at your option) any later version. * + * * + * This program is distributed in the hope that it will be useful, * + * but WITHOUT ANY WARRANTY; without even the implied warranty of * + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * + * GNU General Public License for more details. * + * * + * You should have received a copy of the GNU General Public License * + * along with this program. If not, see . * + *****************************************************************************/ +/*! + * \file + * + * \brief A fluid system with water and a ficitious component, which is to be + * implemented in tutorial exercise 3a, as phases and components. This + * fluid system is to be implemented in exercise 3b. + */ +#ifndef DUMUX_H2O_MYCOMPRESSIBLECOMPONENT_FLUID_SYSTEM_HH +#define DUMUX_H2O_MYCOMPRESSIBLECOMPONENT_FLUID_SYSTEM_HH + +#include +#include + +#include + +// the ficitious component that was created in exercise 3a +#include + +// the binary coefficients corresponding to this fluid system +#include + +namespace Dumux +{ +namespace FluidSystems +{ + +/*! + * \brief A compositional fluid consisting of two liquid phases, + * which are water and a ficitious component from tutorial exercise 3a. + */ +template > > +class H2OMyCompressibleComponent + : public BaseFluidSystem< Scalar, H2OMyCompressibleComponent > +{ + typedef H2OMyCompressibleComponent ThisType; + typedef BaseFluidSystem Base; + +public: + typedef Dumux::MyCompressibleComponent MyCompressibleComponent; + typedef H2OType H2O; + + static const int numPhases = 2; + static const int numComponents = 2; + + static const int wPhaseIdx = 0; // index of the water phase + static const int nPhaseIdx = 1; // index of the NAPL phase + + static const int H2OIdx = 0; + static const int NAPLIdx = 1; + + // export component indices to indicate the main component + // of the corresponding phase at atmospheric pressure 1 bar + // and room temperature 20°C: + static const int wCompIdx = H2OIdx; + static const int nCompIdx = NAPLIdx; + + /*! + * \brief Initialize the fluid system's static parameters generically + * + * If a tabulated H2O component is used, we do our best to create + * tables that always work. + */ + static void init() + { + init(/*tempMin=*/273.15, + /*tempMax=*/623.15, + /*numTemp=*/100, + /*pMin=*/0.0, + /*pMax=*/20e6, + /*numP=*/200); + } + + /*! + * \brief Initialize the fluid system's static parameters using + * problem specific temperature and pressure ranges + * + * \param tempMin The minimum temperature used for tabulation of water [K] + * \param tempMax The maximum temperature used for tabulation of water [K] + * \param nTemp The number of ticks on the temperature axis of the table of water + * \param pressMin The minimum pressure used for tabulation of water [Pa] + * \param pressMax The maximum pressure used for tabulation of water [Pa] + * \param nPress The number of ticks on the pressure axis of the table of water + */ + static void init(Scalar tempMin, Scalar tempMax, unsigned nTemp, + Scalar pressMin, Scalar pressMax, unsigned nPress) + { + if (H2O::isTabulated) { + std::cout << "Initializing tables for the H2O fluid properties (" + << nTemp*nPress + << " entries).\n"; + + H2O::init(tempMin, tempMax, nTemp, + pressMin, pressMax, nPress); + } + } + + + /*! + * \brief Return whether a phase is liquid + * + * \param phaseIdx The index of the fluid phase to consider + */ + static bool isLiquid(int phaseIdx) + { + assert(0 <= phaseIdx && phaseIdx < numPhases); + return true; + } + + static bool isIdealGas(int phaseIdx) + { return H2O::gasIsIdeal() && MyCompressibleComponent::gasIsIdeal(); } + + /*! + * \brief Returns true if and only if a fluid phase is assumed to + * be an ideal mixture. + * + * We define an ideal mixture as a fluid phase where the fugacity + * coefficients of all components times the pressure of the phase + * are indepent on the fluid composition. This assumtion is true + * if Henry's law and Raoult's law apply. If you are unsure what + * this function should return, it is safe to return false. The + * only damage done will be (slightly) increased computation times + * in some cases. + * + * \param phaseIdx The index of the fluid phase to consider + */ + static bool isIdealMixture(int phaseIdx) + { + assert(0 <= phaseIdx && phaseIdx < numPhases); + return true; + } + + /*! + * \brief Returns true if and only if a fluid phase is assumed to + * be compressible. + * + * Compressible means that the partial derivative of the density + * to the fluid pressure is always larger than zero. + * + * \param phaseIdx The index of the fluid phase to consider + */ + static bool isCompressible(int phaseIdx) + { + assert(0 <= phaseIdx && phaseIdx < numPhases); + if (phaseIdx == wPhaseIdx) + // the water component decides for the water phase... + return H2O::liquidIsCompressible(); + + // the NAPL component decides for the napl phase... + return MyCompressibleComponent::liquidIsCompressible(); + } + + /*! + * \brief Return the human readable name of a phase (used in indices) + */ + static std::string phaseName(int phaseIdx) + { + switch (phaseIdx) { + case wPhaseIdx: return "w"; + case nPhaseIdx: return "n"; + }; + DUNE_THROW(Dune::InvalidStateException, "Invalid phase index " << phaseIdx); + } + + /*! + * \brief Return the human readable name of a component (used in indices) + */ + static std::string componentName(int compIdx) + { + switch (compIdx) { + case H2OIdx: return H2O::name(); + case NAPLIdx: return MyCompressibleComponent::name(); + }; + DUNE_THROW(Dune::InvalidStateException, "Invalid component index " << compIdx); + } + + /*! + * \brief Return the molar mass of a component in [kg/mol]. + */ + static Scalar molarMass(int compIdx) + { + switch (compIdx) { + case H2OIdx: return H2O::molarMass(); + case NAPLIdx: return MyCompressibleComponent::molarMass(); + }; + DUNE_THROW(Dune::InvalidStateException, "Invalid component index " << compIdx); + } + + /*! + * \brief Given all mole fractions in a phase, return the phase + * density [kg/m^3]. + */ + using Base::density; + template + static Scalar density(const FluidState &fluidState, int phaseIdx) + { + assert(0 <= phaseIdx && phaseIdx < numPhases); + if (phaseIdx == wPhaseIdx) { + // See: doctoral thesis of Steffen Ochs 2007 + // Steam injection into saturated porous media : process analysis including experimental and numerical investigations + // http://elib.uni-stuttgart.de/bitstream/11682/271/1/Diss_Ochs_OPUS.pdf + Scalar rholH2O = H2O::liquidDensity(fluidState.temperature(phaseIdx), fluidState.pressure(phaseIdx)); + Scalar clH2O = rholH2O/H2O::molarMass(); + Scalar x_H2O = fluidState.moleFraction(wPhaseIdx, H2OIdx); + Scalar x_myComp = fluidState.moleFraction(wPhaseIdx, NAPLIdx); + + // return composition-dependent water phase density + return clH2O*(H2O::molarMass()*x_H2O + MyCompressibleComponent::molarMass()*x_myComp); + } + else { + // assume the density of the fictious component to be independent of the composition + Scalar pressure = MyCompressibleComponent::liquidIsCompressible()?fluidState.pressure(phaseIdx):1e100; + return MyCompressibleComponent::liquidDensity(fluidState.temperature(phaseIdx), pressure); + } + } + + /*! + * \brief Return the viscosity of a phase. + */ + using Base::viscosity; + template + static Scalar viscosity(const FluidState &fluidState, + int phaseIdx) + { + assert(0 <= phaseIdx && phaseIdx < numPhases); + if (phaseIdx == wPhaseIdx) { + // assume pure water viscosity + return H2O::liquidViscosity(fluidState.temperature(phaseIdx), + fluidState.pressure(phaseIdx)); + } + else { + // assume pure NAPL viscosity + return MyCompressibleComponent::liquidViscosity(fluidState.temperature(phaseIdx), fluidState.pressure(phaseIdx)); + } + } + + using Base::diffusionCoefficient; + template + static Scalar diffusionCoefficient(const FluidState &fluidState, int phaseIdx, int compIdx) + { + DUNE_THROW(Dune::NotImplemented, "Diffusion coefficients"); + } + + /*! + * \brief Given a phase's composition, temperature and pressure, + * return the binary diffusion coefficient \f$\mathrm{[m^2/s]}\f$ for components + * \f$\mathrm{i}\f$ and \f$\mathrm{j}\f$ in this phase. + * \param fluidState The fluid state + * \param paramCache mutable parameters + * \param phaseIdx Index of the fluid phase + * \param compIIdx Index of the component i + * \param compJIdx Index of the component j + */ + using Base::binaryDiffusionCoefficient; + template + static Scalar binaryDiffusionCoefficient(const FluidState &fluidState, + int phaseIdx, + int compIIdx, + int compJIdx) + { + assert(0 <= phaseIdx && phaseIdx < numPhases); + assert(0 <= compIdx && compIdx < numComponents); + + const Scalar T = fluidState.temperature(phaseIdx); + const Scalar p = fluidState.pressure(phaseIdx); + + // we assume the diffusion coefficient to be the same in both phases + return Dumux::BinaryCoeff::H2O_MyCompressibleComponent::liquidDiffCoeff(T, p); + } + + /* Henry coefficients + */ + template + static Scalar henryCoefficient(const FluidState &fluidState, + int phaseIdx, + int compIdx) + { + assert(0 <= phaseIdx && phaseIdx < numPhases); + assert(0 <= compIdx && compIdx < numComponents); + + const Scalar T = fluidState.temperature(phaseIdx); + const Scalar p = fluidState.pressure(phaseIdx); + + if (compIdx == NAPLIdx && phaseIdx == wPhaseIdx) + return Dumux::BinaryCoeff::H2O_MyCompressibleComponent::henryMyCompressibleComponentInWater(T)/p; + + else if (phaseIdx == nPhaseIdx && compIdx == H2OIdx) + return Dumux::BinaryCoeff::H2O_MyCompressibleComponent::henryWaterInMyCompressibleComponent(T)/p; + + else + DUNE_THROW(Dune::InvalidStateException, "non-existent henry coefficient for phase index " << phaseIdx + << " and component index " << compIdx); + } + + using Base::fugacityCoefficient; + /*! + * \brief Returns the fugacity coefficient \f$\mathrm{[-]}\f$ of a component in a + * phase. + * + * In this case, things are actually pretty simple. We have an ideal + * solution. Thus, the fugacity coefficient is 1 in the gas phase + * (fugacity equals the partial pressure of the component in the gas phase + * respectively in the liquid phases it is the inverse of the + * Henry coefficients scaled by pressure + * \param fluidState The fluid state + * \param phaseIdx The index of the phase + * \param compIdx The index of the component + */ + template + static Scalar fugacityCoefficient(const FluidState &fluidState, + int phaseIdx, + int compIdx) + { + assert(0 <= phaseIdx && phaseIdx < numPhases); + assert(0 <= compIdx && compIdx < numComponents); + + Scalar T = fluidState.temperature(phaseIdx); + Scalar p = fluidState.pressure(phaseIdx); + + if (phaseIdx == wPhaseIdx) { + if (compIdx == H2OIdx) + return H2O::vaporPressure(T)/p; + else if (compIdx == NAPLIdx) + return Dumux::BinaryCoeff::H2O_MyCompressibleComponent::henryMyCompressibleComponentInWater(T)/p; + } + + // for the NAPL phase, we assume currently that nothing is + // dissolved. this means that the affinity of the NAPL + // component to the NAPL phase is much higher than for the + // other components, i.e. the fugacity coefficient is much + // smaller. + Scalar phiNapl = MyCompressibleComponent::vaporPressure(T)/p; + if (compIdx == NAPLIdx) + return phiNapl; + else + return 1e6*phiNapl; + } + + template + static Scalar kelvinVaporPressure(const FluidState &fluidState, + const int phaseIdx, + const int compIdx) + { + DUNE_THROW(Dune::NotImplemented, "FluidSystems::H2OMyCompressibleComponent::kelvinVaporPressure()"); + } + + /* partial pressures in the gas phase, taken from saturation vapor pressures + */ + template + static Scalar partialPressureGas(const FluidState &fluidState, int phaseIdx, + int compIdx) + { + assert(0 <= compIdx && compIdx < numComponents); + + const Scalar T = fluidState.temperature(phaseIdx); + if (compIdx == NAPLIdx) + return MyCompressibleComponent::vaporPressure(T); + else if (compIdx == H2OIdx) + return H2O::vaporPressure(T); + else + DUNE_THROW(Dune::InvalidStateException, "non-existent component index " << compIdx); + } + + /* inverse vapor pressures, taken from inverse saturation vapor pressures + */ + template + static Scalar inverseVaporPressureCurve(const FluidState &fluidState, + int phaseIdx, + int compIdx) + { + assert(0 <= compIdx && compIdx < numComponents); + + const Scalar pressure = fluidState.pressure(phaseIdx); + if (compIdx == NAPLIdx) + return MyCompressibleComponent::vaporTemperature(pressure); + else if (compIdx == H2OIdx) + return H2O::vaporTemperature(pressure); + else + DUNE_THROW(Dune::InvalidStateException, "non-existent component index " << compIdx); + } + + + + /*! + * \brief Given all mole fractions in a phase, return the specific + * phase enthalpy [J/kg]. + */ + using Base::enthalpy; + template + static Scalar enthalpy(const FluidState &fluidState, + int phaseIdx) + { + assert(0 <= phaseIdx && phaseIdx < numPhases); + if (phaseIdx == wPhaseIdx) { + return H2O::liquidEnthalpy(fluidState.temperature(phaseIdx), fluidState.pressure(phaseIdx)); + } + else { + return MyCompressibleComponent::liquidEnthalpy(fluidState.temperature(phaseIdx), fluidState.pressure(phaseIdx)); + } + DUNE_THROW(Dune::InvalidStateException, "Invalid phase index " << phaseIdx); + } + + using Base::heatCapacity; + template + static Scalar heatCapacity(const FluidState &fluidState, + int phaseIdx) + { + DUNE_THROW(Dune::NotImplemented, "FluidSystems::H2ONAPL::heatCapacity()"); + } + + using Base::thermalConductivity; + template + static Scalar thermalConductivity(const FluidState &fluidState, + int phaseIdx) + { + assert(0 <= phaseIdx && phaseIdx < numPhases); + + const Scalar temperature = fluidState.temperature(phaseIdx) ; + const Scalar pressure = fluidState.pressure(phaseIdx); + if (phaseIdx == wPhaseIdx) + { + return H2O::liquidThermalConductivity(temperature, pressure); + } + else + { + return MyCompressibleComponent::liquidThermalConductivity(temperature, pressure); + } + DUNE_THROW(Dune::InvalidStateException, "Invalid phase index " << phaseIdx); + } + +private: + +}; +} // end namespace FluidSystems +} // end namespace Dumux + +#endif diff --git a/tutorial/solution/ex3/mycompressiblecomponent.hh b/tutorial/solution/ex3/mycompressiblecomponent.hh new file mode 100644 index 0000000000000000000000000000000000000000..6f349b541d32dd890585546f5da7a906bc910e78 --- /dev/null +++ b/tutorial/solution/ex3/mycompressiblecomponent.hh @@ -0,0 +1,108 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/***************************************************************************** + * See the file COPYING for full copying permissions. * + * * + * This program is free software: you can redistribute it and/or modify * + * it under the terms of the GNU General Public License as published by * + * the Free Software Foundation, either version 2 of the License, or * + * (at your option) any later version. * + * * + * This program is distributed in the hope that it will be useful, * + * but WITHOUT ANY WARRANTY; without even the implied warranty of * + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * + * GNU General Public License for more details. * + * * + * You should have received a copy of the GNU General Public License * + * along with this program. If not, see . * + *****************************************************************************/ +/*! + * \file + * \ingroup Components + * \brief A ficitious component to be implemented in exercise 3. + */ +#ifndef DUMUX_MYCOMPRESSIBLECOMPONENT_HH +#define DUMUX_MYCOMPRESSIBLECOMPONENT_HH + +#include +#include + + +namespace Dumux +{ +/*! + * \ingroup Components + * \brief A ficitious component to be implemented in exercise 3. + * + * \tparam Scalar The type used for scalar values + */ +template +class MyCompressibleComponent : public Component > +{ +public: + /*! + * \brief A human readable name for MyCompressibleComponent. + */ + static std::string name() + { return "MyCompressibleComponent"; } + + /*! + * \brief The molar mass in \f$\mathrm{[kg/mol]}\f$ of MyCompressibleComponent. + */ + static Scalar molarMass() + { + return 131.39e-3; // [kg/mol] + } + + /*! + * \brief The density of MyCompressibleComponent at a given pressure and temperature \f$\mathrm{[kg/m^3]}\f$. + * + * \param temperature temperature of component in \f$\mathrm{[K]}\f$ + * \param pressure pressure of component in \f$\mathrm{[Pa]}\f$ + */ + static Scalar liquidDensity(Scalar temperature, Scalar pressure) + { + static const Scalar rho_min = 1440; + static const Scalar rho_max = 1480; + static const Scalar k = 5e-7; + + using std::exp; + return rho_min + (rho_max - rho_min)/(1 + rho_min*exp(-1.0*k*(rho_max - rho_min)*pressure)); // [kg/m^3] + } + + /*! + * \brief The dynamic viscosity \f$\mathrm{[Pa*s]}\f$ of MyCompressibleComponent. + * + * \param temperature temperature of component in \f$\mathrm{[K]}\f$ + * \param pressure pressure of component in \f$\mathrm{[Pa]}\f$ + */ + static Scalar liquidViscosity(Scalar temperature, Scalar pressure) + { + return 5.7e-4;// [Pa*s] + } + + /***************************************************************** + * The function below is implemented in the scope of exercise 3b * + *****************************************************************/ + + /*! + * \brief The vapor pressure in \f$\mathrm{[Pa]}\f$ of MyCompressibleComponent + * at a given temperature. + * + * \param T temperature of component in \f$\mathrm{[K]}\f$ + */ + static Scalar vaporPressure(Scalar T) + { + return 3900; // [Pa] (at 20C) + } + + /*! + * \brief Returns true if the liquid phase is assumed to be compressible + */ + static bool liquidIsCompressible() + { return false; } +}; + +} // end namespace + +#endif diff --git a/tutorial/solution/ex3/myincompressiblecomponent.hh b/tutorial/solution/ex3/myincompressiblecomponent.hh new file mode 100644 index 0000000000000000000000000000000000000000..f8af80ea8a5da8c2cce1fe3f7e37fc2f0bc56129 --- /dev/null +++ b/tutorial/solution/ex3/myincompressiblecomponent.hh @@ -0,0 +1,82 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/***************************************************************************** + * See the file COPYING for full copying permissions. * + * * + * This program is free software: you can redistribute it and/or modify * + * it under the terms of the GNU General Public License as published by * + * the Free Software Foundation, either version 2 of the License, or * + * (at your option) any later version. * + * * + * This program is distributed in the hope that it will be useful, * + * but WITHOUT ANY WARRANTY; without even the implied warranty of * + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * + * GNU General Public License for more details. * + * * + * You should have received a copy of the GNU General Public License * + * along with this program. If not, see . * + *****************************************************************************/ +/*! + * \file + * \ingroup Components + * \brief A ficitious component to be implemented in tutorial exercise 3. + */ +#ifndef DUMUX_MYINCOMPRESSIBLECOMPONENT_HH +#define DUMUX_MYINCOMPRESSIBLECOMPONENT_HH + +#include +#include + + +namespace Dumux +{ +/*! + * \ingroup Components + * \brief A ficitious component to be implemented in exercise 3. + * + * \tparam Scalar The type used for scalar values + */ +template +class MyIncompressibleComponent : public Component > +{ +public: + /*! + * \brief A human readable name for MyIncompressibleComponent. + */ + static std::string name() + { return "MyIncompressibleComponent"; } + + /*! + * \brief The molar mass in \f$\mathrm{[kg/mol]}\f$ of MyIncompressibleComponent. + */ + static Scalar molarMass() + { + return 131.39e-3; // [kg/mol] + } + + /*! + * \brief The density of MyIncompressibleComponent at a given pressure and temperature \f$\mathrm{[kg/m^3]}\f$. + * + * \param temperature temperature of component in \f$\mathrm{[K]}\f$ + * \param pressure pressure of component in \f$\mathrm{[Pa]}\f$ + */ + static Scalar liquidDensity(Scalar temperature, Scalar pressure) + { + return 1460.0; // [kg/m^3] + } + + /*! + * \brief The dynamic viscosity \f$\mathrm{[Pa*s]}\f$ of MyIncompressibleComponent. + * + * \param temperature temperature of component in \f$\mathrm{[K]}\f$ + * \param pressure pressure of component in \f$\mathrm{[Pa]}\f$ + */ + static Scalar liquidViscosity(Scalar temperature, Scalar pressure) + { + return 5.7e-4;// [Pa*s] + } +}; + +} // end namespace + +#endif diff --git a/tutorial/tutorial_implicit/CMakeLists.txt b/tutorial/tutorial_implicit/CMakeLists.txt new file mode 100644 index 0000000000000000000000000000000000000000..e73f579367bbf568540c7475182a06345f456922 --- /dev/null +++ b/tutorial/tutorial_implicit/CMakeLists.txt @@ -0,0 +1,11 @@ +add_input_file_links() + +add_dumux_test(tutorial_implicit tutorial_implicit tutorial_implicit.cc + ${CMAKE_CURRENT_BINARY_DIR}/tutorial_implicit) + +#install sources +install(FILES +tutorial_implicit.cc +tutorialproblem_implicit.hh +tutorialspatialparams_implicit.hh +DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dumux/tutorial/tutorial_implicit) diff --git a/tutorial/solutions_implicit/ex1b_tutorial_implicit.input.diff b/tutorial/tutorial_implicit/solutions_implicit/ex1b_tutorial_implicit.input.diff similarity index 100% rename from tutorial/solutions_implicit/ex1b_tutorial_implicit.input.diff rename to tutorial/tutorial_implicit/solutions_implicit/ex1b_tutorial_implicit.input.diff diff --git a/tutorial/solutions_implicit/ex1b_tutorialproblem_implicit.diff b/tutorial/tutorial_implicit/solutions_implicit/ex1b_tutorialproblem_implicit.diff similarity index 100% rename from tutorial/solutions_implicit/ex1b_tutorialproblem_implicit.diff rename to tutorial/tutorial_implicit/solutions_implicit/ex1b_tutorialproblem_implicit.diff diff --git a/tutorial/solutions_implicit/ex1c_tutorialproblem_implicit.diff b/tutorial/tutorial_implicit/solutions_implicit/ex1c_tutorialproblem_implicit.diff similarity index 100% rename from tutorial/solutions_implicit/ex1c_tutorialproblem_implicit.diff rename to tutorial/tutorial_implicit/solutions_implicit/ex1c_tutorialproblem_implicit.diff diff --git a/tutorial/solutions_implicit/ex1d_tutorialproblem_implicit.diff b/tutorial/tutorial_implicit/solutions_implicit/ex1d_tutorialproblem_implicit.diff similarity index 100% rename from tutorial/solutions_implicit/ex1d_tutorialproblem_implicit.diff rename to tutorial/tutorial_implicit/solutions_implicit/ex1d_tutorialproblem_implicit.diff diff --git a/tutorial/solutions_implicit/ex1e_tutorialproblem_implicit.diff b/tutorial/tutorial_implicit/solutions_implicit/ex1e_tutorialproblem_implicit.diff similarity index 100% rename from tutorial/solutions_implicit/ex1e_tutorialproblem_implicit.diff rename to tutorial/tutorial_implicit/solutions_implicit/ex1e_tutorialproblem_implicit.diff diff --git a/tutorial/solutions_implicit/ex1f_tutorialspatialparams_implicit.diff b/tutorial/tutorial_implicit/solutions_implicit/ex1f_tutorialspatialparams_implicit.diff similarity index 100% rename from tutorial/solutions_implicit/ex1f_tutorialspatialparams_implicit.diff rename to tutorial/tutorial_implicit/solutions_implicit/ex1f_tutorialspatialparams_implicit.diff diff --git a/tutorial/solutions_implicit/ex1g_tutorialspatialparams_implicit.diff b/tutorial/tutorial_implicit/solutions_implicit/ex1g_tutorialspatialparams_implicit.diff similarity index 100% rename from tutorial/solutions_implicit/ex1g_tutorialspatialparams_implicit.diff rename to tutorial/tutorial_implicit/solutions_implicit/ex1g_tutorialspatialparams_implicit.diff diff --git a/tutorial/solutions_implicit/ex2_tutorial_implicit_input.diff b/tutorial/tutorial_implicit/solutions_implicit/ex2_tutorial_implicit_input.diff similarity index 100% rename from tutorial/solutions_implicit/ex2_tutorial_implicit_input.diff rename to tutorial/tutorial_implicit/solutions_implicit/ex2_tutorial_implicit_input.diff diff --git a/tutorial/solutions_implicit/ex2_tutorialproblem_implicit.hh b/tutorial/tutorial_implicit/solutions_implicit/ex2_tutorialproblem_implicit.hh similarity index 100% rename from tutorial/solutions_implicit/ex2_tutorialproblem_implicit.hh rename to tutorial/tutorial_implicit/solutions_implicit/ex2_tutorialproblem_implicit.hh diff --git a/tutorial/solutions_implicit/ex2_tutorialspatialparams_implicit.hh b/tutorial/tutorial_implicit/solutions_implicit/ex2_tutorialspatialparams_implicit.hh similarity index 100% rename from tutorial/solutions_implicit/ex2_tutorialspatialparams_implicit.hh rename to tutorial/tutorial_implicit/solutions_implicit/ex2_tutorialspatialparams_implicit.hh diff --git a/tutorial/solutions_implicit/ex3_tutorial_implicit.input b/tutorial/tutorial_implicit/solutions_implicit/ex3_tutorial_implicit.input similarity index 100% rename from tutorial/solutions_implicit/ex3_tutorial_implicit.input rename to tutorial/tutorial_implicit/solutions_implicit/ex3_tutorial_implicit.input diff --git a/tutorial/solutions_implicit/ex3_tutorialproblem_implicit.diff b/tutorial/tutorial_implicit/solutions_implicit/ex3_tutorialproblem_implicit.diff similarity index 100% rename from tutorial/solutions_implicit/ex3_tutorialproblem_implicit.diff rename to tutorial/tutorial_implicit/solutions_implicit/ex3_tutorialproblem_implicit.diff diff --git a/tutorial/solutions_implicit/ex3_tutorialspatialparams_implicit.diff b/tutorial/tutorial_implicit/solutions_implicit/ex3_tutorialspatialparams_implicit.diff similarity index 100% rename from tutorial/solutions_implicit/ex3_tutorialspatialparams_implicit.diff rename to tutorial/tutorial_implicit/solutions_implicit/ex3_tutorialspatialparams_implicit.diff diff --git a/tutorial/solutions_implicit/ex4_benzene.hh b/tutorial/tutorial_implicit/solutions_implicit/ex4_benzene.hh similarity index 100% rename from tutorial/solutions_implicit/ex4_benzene.hh rename to tutorial/tutorial_implicit/solutions_implicit/ex4_benzene.hh diff --git a/tutorial/solutions_implicit/ex4_tutorialproblem_implicit.diff b/tutorial/tutorial_implicit/solutions_implicit/ex4_tutorialproblem_implicit.diff similarity index 100% rename from tutorial/solutions_implicit/ex4_tutorialproblem_implicit.diff rename to tutorial/tutorial_implicit/solutions_implicit/ex4_tutorialproblem_implicit.diff diff --git a/tutorial/solutions_implicit/ex5_tutorial_implicit.input b/tutorial/tutorial_implicit/solutions_implicit/ex5_tutorial_implicit.input similarity index 100% rename from tutorial/solutions_implicit/ex5_tutorial_implicit.input rename to tutorial/tutorial_implicit/solutions_implicit/ex5_tutorial_implicit.input diff --git a/tutorial/solutions_implicit/ex5_tutorialproblem_implicit.hh b/tutorial/tutorial_implicit/solutions_implicit/ex5_tutorialproblem_implicit.hh similarity index 100% rename from tutorial/solutions_implicit/ex5_tutorialproblem_implicit.hh rename to tutorial/tutorial_implicit/solutions_implicit/ex5_tutorialproblem_implicit.hh diff --git a/tutorial/solutions_implicit/ex5_tutorialspatialparams_implicit.hh b/tutorial/tutorial_implicit/solutions_implicit/ex5_tutorialspatialparams_implicit.hh similarity index 100% rename from tutorial/solutions_implicit/ex5_tutorialspatialparams_implicit.hh rename to tutorial/tutorial_implicit/solutions_implicit/ex5_tutorialspatialparams_implicit.hh diff --git a/tutorial/tutorial_implicit.cc b/tutorial/tutorial_implicit/tutorial_implicit.cc similarity index 100% rename from tutorial/tutorial_implicit.cc rename to tutorial/tutorial_implicit/tutorial_implicit.cc diff --git a/tutorial/tutorial_implicit.input b/tutorial/tutorial_implicit/tutorial_implicit.input similarity index 100% rename from tutorial/tutorial_implicit.input rename to tutorial/tutorial_implicit/tutorial_implicit.input diff --git a/tutorial/tutorialproblem_implicit.hh b/tutorial/tutorial_implicit/tutorialproblem_implicit.hh similarity index 100% rename from tutorial/tutorialproblem_implicit.hh rename to tutorial/tutorial_implicit/tutorialproblem_implicit.hh diff --git a/tutorial/tutorialspatialparams_implicit.hh b/tutorial/tutorial_implicit/tutorialspatialparams_implicit.hh similarity index 100% rename from tutorial/tutorialspatialparams_implicit.hh rename to tutorial/tutorial_implicit/tutorialspatialparams_implicit.hh diff --git a/tutorial/tutorial_sequential/CMakeLists.txt b/tutorial/tutorial_sequential/CMakeLists.txt new file mode 100644 index 0000000000000000000000000000000000000000..d4ddb99d63d52a1527436b5cb6636e0b8db157d5 --- /dev/null +++ b/tutorial/tutorial_sequential/CMakeLists.txt @@ -0,0 +1,11 @@ +add_input_file_links() + +add_dumux_test(tutorial_sequential tutorial_sequential tutorial_sequential.cc + ${CMAKE_CURRENT_BINARY_DIR}/tutorial_sequential) + +#install sources +install(FILES +tutorial_sequential.cc +tutorialproblem_sequential.hh +tutorialspatialparams_sequential.hh +DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dumux/tutorial/tutorial_sequential) diff --git a/tutorial/solutions_sequential/Ex1_a_input.diff b/tutorial/tutorial_sequential/solutions_sequential/Ex1_a_input.diff similarity index 100% rename from tutorial/solutions_sequential/Ex1_a_input.diff rename to tutorial/tutorial_sequential/solutions_sequential/Ex1_a_input.diff diff --git a/tutorial/solutions_sequential/Ex1_a_problem.diff b/tutorial/tutorial_sequential/solutions_sequential/Ex1_a_problem.diff similarity index 100% rename from tutorial/solutions_sequential/Ex1_a_problem.diff rename to tutorial/tutorial_sequential/solutions_sequential/Ex1_a_problem.diff diff --git a/tutorial/solutions_sequential/Ex1_b_input.diff b/tutorial/tutorial_sequential/solutions_sequential/Ex1_b_input.diff similarity index 100% rename from tutorial/solutions_sequential/Ex1_b_input.diff rename to tutorial/tutorial_sequential/solutions_sequential/Ex1_b_input.diff diff --git a/tutorial/solutions_sequential/Ex1_b_problem.diff b/tutorial/tutorial_sequential/solutions_sequential/Ex1_b_problem.diff similarity index 100% rename from tutorial/solutions_sequential/Ex1_b_problem.diff rename to tutorial/tutorial_sequential/solutions_sequential/Ex1_b_problem.diff diff --git a/tutorial/solutions_sequential/Ex1_c_problem.diff b/tutorial/tutorial_sequential/solutions_sequential/Ex1_c_problem.diff similarity index 100% rename from tutorial/solutions_sequential/Ex1_c_problem.diff rename to tutorial/tutorial_sequential/solutions_sequential/Ex1_c_problem.diff diff --git a/tutorial/solutions_sequential/Ex1_d_input.diff b/tutorial/tutorial_sequential/solutions_sequential/Ex1_d_input.diff similarity index 100% rename from tutorial/solutions_sequential/Ex1_d_input.diff rename to tutorial/tutorial_sequential/solutions_sequential/Ex1_d_input.diff diff --git a/tutorial/solutions_sequential/Ex1_d_problem.diff b/tutorial/tutorial_sequential/solutions_sequential/Ex1_d_problem.diff similarity index 100% rename from tutorial/solutions_sequential/Ex1_d_problem.diff rename to tutorial/tutorial_sequential/solutions_sequential/Ex1_d_problem.diff diff --git a/tutorial/solutions_sequential/Ex1_e_input.diff b/tutorial/tutorial_sequential/solutions_sequential/Ex1_e_input.diff similarity index 100% rename from tutorial/solutions_sequential/Ex1_e_input.diff rename to tutorial/tutorial_sequential/solutions_sequential/Ex1_e_input.diff diff --git a/tutorial/solutions_sequential/Ex1_e_spatialparams.diff b/tutorial/tutorial_sequential/solutions_sequential/Ex1_e_spatialparams.diff similarity index 100% rename from tutorial/solutions_sequential/Ex1_e_spatialparams.diff rename to tutorial/tutorial_sequential/solutions_sequential/Ex1_e_spatialparams.diff diff --git a/tutorial/solutions_sequential/Ex2_tutorial_sequential.input.diff b/tutorial/tutorial_sequential/solutions_sequential/Ex2_tutorial_sequential.input.diff similarity index 100% rename from tutorial/solutions_sequential/Ex2_tutorial_sequential.input.diff rename to tutorial/tutorial_sequential/solutions_sequential/Ex2_tutorial_sequential.input.diff diff --git a/tutorial/solutions_sequential/Ex3tutorial_sequential.input.diff b/tutorial/tutorial_sequential/solutions_sequential/Ex3tutorial_sequential.input.diff similarity index 100% rename from tutorial/solutions_sequential/Ex3tutorial_sequential.input.diff rename to tutorial/tutorial_sequential/solutions_sequential/Ex3tutorial_sequential.input.diff diff --git a/tutorial/solutions_sequential/Ex3tutorialproblem_sequential.diff b/tutorial/tutorial_sequential/solutions_sequential/Ex3tutorialproblem_sequential.diff similarity index 100% rename from tutorial/solutions_sequential/Ex3tutorialproblem_sequential.diff rename to tutorial/tutorial_sequential/solutions_sequential/Ex3tutorialproblem_sequential.diff diff --git a/tutorial/solutions_sequential/Ex3tutorialspatialparams_sequential.diff b/tutorial/tutorial_sequential/solutions_sequential/Ex3tutorialspatialparams_sequential.diff similarity index 100% rename from tutorial/solutions_sequential/Ex3tutorialspatialparams_sequential.diff rename to tutorial/tutorial_sequential/solutions_sequential/Ex3tutorialspatialparams_sequential.diff diff --git a/tutorial/solutions_sequential/ex2tutorialproblem_sequential.hh b/tutorial/tutorial_sequential/solutions_sequential/ex2tutorialproblem_sequential.hh similarity index 100% rename from tutorial/solutions_sequential/ex2tutorialproblem_sequential.hh rename to tutorial/tutorial_sequential/solutions_sequential/ex2tutorialproblem_sequential.hh diff --git a/tutorial/solutions_sequential/ex2tutorialspatialparams_sequential.hh b/tutorial/tutorial_sequential/solutions_sequential/ex2tutorialspatialparams_sequential.hh similarity index 100% rename from tutorial/solutions_sequential/ex2tutorialspatialparams_sequential.hh rename to tutorial/tutorial_sequential/solutions_sequential/ex2tutorialspatialparams_sequential.hh diff --git a/tutorial/solutions_sequential/ex4_benzene.hh b/tutorial/tutorial_sequential/solutions_sequential/ex4_benzene.hh similarity index 100% rename from tutorial/solutions_sequential/ex4_benzene.hh rename to tutorial/tutorial_sequential/solutions_sequential/ex4_benzene.hh diff --git a/tutorial/solutions_sequential/ex5_tutorial_sequential.input b/tutorial/tutorial_sequential/solutions_sequential/ex5_tutorial_sequential.input similarity index 100% rename from tutorial/solutions_sequential/ex5_tutorial_sequential.input rename to tutorial/tutorial_sequential/solutions_sequential/ex5_tutorial_sequential.input diff --git a/tutorial/solutions_sequential/ex5_tutorialproblem_sequential.hh b/tutorial/tutorial_sequential/solutions_sequential/ex5_tutorialproblem_sequential.hh similarity index 100% rename from tutorial/solutions_sequential/ex5_tutorialproblem_sequential.hh rename to tutorial/tutorial_sequential/solutions_sequential/ex5_tutorialproblem_sequential.hh diff --git a/tutorial/solutions_sequential/ex5_tutorialspatialparams_sequential.hh b/tutorial/tutorial_sequential/solutions_sequential/ex5_tutorialspatialparams_sequential.hh similarity index 100% rename from tutorial/solutions_sequential/ex5_tutorialspatialparams_sequential.hh rename to tutorial/tutorial_sequential/solutions_sequential/ex5_tutorialspatialparams_sequential.hh diff --git a/tutorial/tutorial_sequential.cc b/tutorial/tutorial_sequential/tutorial_sequential.cc similarity index 100% rename from tutorial/tutorial_sequential.cc rename to tutorial/tutorial_sequential/tutorial_sequential.cc diff --git a/tutorial/tutorial_sequential.input b/tutorial/tutorial_sequential/tutorial_sequential.input similarity index 100% rename from tutorial/tutorial_sequential.input rename to tutorial/tutorial_sequential/tutorial_sequential.input diff --git a/tutorial/tutorialproblem_sequential.hh b/tutorial/tutorial_sequential/tutorialproblem_sequential.hh similarity index 100% rename from tutorial/tutorialproblem_sequential.hh rename to tutorial/tutorial_sequential/tutorialproblem_sequential.hh diff --git a/tutorial/tutorialspatialparams_sequential.hh b/tutorial/tutorial_sequential/tutorialspatialparams_sequential.hh similarity index 100% rename from tutorial/tutorialspatialparams_sequential.hh rename to tutorial/tutorial_sequential/tutorialspatialparams_sequential.hh