diff --git a/exercises/exercise-fluidsystem/2p2cproblem.hh b/exercises/exercise-fluidsystem/2p2cproblem.hh
index 13eb748b8a9f18e116d75d9b8746fed6a6ecffe9..8cff7cc77480322f0a88bff3ba194ea3c8a9313d 100644
--- a/exercises/exercise-fluidsystem/2p2cproblem.hh
+++ b/exercises/exercise-fluidsystem/2p2cproblem.hh
@@ -138,7 +138,6 @@ public:
* \brief Specifies which kind of boundary condition should be
* used for which equation on a given boundary segment.
*
- * \param bcTypes The boundary types for the conservation equations
* \param globalPos The position for which the bc type should be evaluated
*/
BoundaryTypes boundaryTypesAtPos(const GlobalPosition &globalPos) const
@@ -169,12 +168,7 @@ public:
* \brief Evaluate 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 position of the integration point of the boundary segment.
- *
- * For this method, the \a values parameter stores the mass flux
- * in normal direction of each phase. Negative values mean influx.
*/
PrimaryVariables neumannAtPos(const GlobalPosition &globalPos) const
{
@@ -208,7 +202,6 @@ public:
/*!
* \brief Evaluate the initial value for a control volume.
- *
* \param globalPos The position for which the initial condition should be evaluated
*
* For this method, the \a values parameter stores primary
@@ -235,9 +228,6 @@ public:
/*!
* \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
*/
PrimaryVariables sourceAtPos(const GlobalPosition &globalPos) const
diff --git a/exercises/exercise-fluidsystem/2pproblem.hh b/exercises/exercise-fluidsystem/2pproblem.hh
index dc7347aec91ef1b3c109034ef085357e1fa2f89e..c5decbfe9a18b8bde3565693ea6ba9c98d40629a 100644
--- a/exercises/exercise-fluidsystem/2pproblem.hh
+++ b/exercises/exercise-fluidsystem/2pproblem.hh
@@ -59,6 +59,9 @@
// The two-phase immiscible fluid system
#include <dumux/material/fluidsystems/2pimmiscible.hh>
+// The interface to create plots during simulation
+#include <dumux/io/gnuplotinterface.hh>
+
namespace Dumux{
// Forward declaration of the problem class
template <class TypeTag> class ExerciseFluidsystemProblemTwoP;
@@ -129,6 +132,7 @@ class ExerciseFluidsystemProblemTwoP : public PorousMediumFlowProblem<TypeTag>
using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry)::LocalView;
using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
+ using FluidState = typename GET_PROP_TYPE(TypeTag, FluidState);
enum {
waterPressureIdx = Indices::pressureIdx,
@@ -156,6 +160,10 @@ public:
// set the depth of the bottom of the reservoir
depthBOR_ = this->fvGridGeometry().bBoxMax()[dimWorld-1];
+
+ // plot density over pressure of the phase consisting of your component
+ if(getParam<bool>("Output.PlotDensity"))
+ plotDensity_();
}
/*!
@@ -178,7 +186,6 @@ public:
* \brief Specifies which kind of boundary condition should be
* used for which equation on a given boundary segment.
*
- * \param bcTypes The boundary types for the conservation equations
* \param globalPos The position for which the bc type should be evaluated
*/
BoundaryTypes boundaryTypesAtPos(const GlobalPosition &globalPos) const
@@ -209,8 +216,6 @@ public:
* \brief Evaluate 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 position of the integration point of the boundary segment.
*
* For this method, the \a values parameter stores the mass flux
@@ -270,8 +275,6 @@ public:
/*!
* \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
*/
PrimaryVariables sourceAtPos(const GlobalPosition& globalPos) const
@@ -282,8 +285,37 @@ public:
}
private:
+ void plotDensity_()
+ {
+ FluidState fluidState;
+ fluidState.setTemperature(temperature());
+ int numberOfPoints = 100;
+ Scalar xMin = 1e4;
+ Scalar xMax = 1e7;
+ Scalar spacing = std::pow((xMax/xMin), 1.0/(numberOfPoints-1));
+ std::vector<double> x(numberOfPoints);
+ std::vector<double> y(numberOfPoints);
+ for (int i=0; i<numberOfPoints; ++i)
+ x[i] = xMin*std::pow(spacing,i);
+ for (int i=0; i<numberOfPoints; ++i)
+ {
+ fluidState.setPressure(FluidSystem::phase1Idx, x[i]);
+ y[i] = FluidSystem::density(fluidState, FluidSystem::phase1Idx);
+ }
+ gnuplot_.resetPlot();
+ gnuplot_.setXRange(xMin, xMax);
+ gnuplot_.setOption("set logscale x 10");
+ gnuplot_.setOption("set key left top");
+ gnuplot_.setYRange(1440, 1480);
+ gnuplot_.setXlabel("pressure [Pa]");
+ gnuplot_.setYlabel("density [kg/m^3]");
+ gnuplot_.addDataSetToPlot(x, y, "YourComponentPhase_density.dat", "w l t 'Phase consisting of your component'");
+ gnuplot_.plot("YourComponentPhase_density");
+ }
+
Scalar eps_; //! small epsilon value
Scalar depthBOR_; //! depth at the bottom of the reservoir
+ Dumux::GnuplotInterface<double> gnuplot_; //! collects data for plotting
};
}
diff --git a/exercises/exercise-fluidsystem/README.md b/exercises/exercise-fluidsystem/README.md
index 75694288bc5170484a2a981e52059cf5929d976a..c9db143d4e66746da540bfb0cb088c3fe75e61aa 100644
--- a/exercises/exercise-fluidsystem/README.md
+++ b/exercises/exercise-fluidsystem/README.md
@@ -1,6 +1,6 @@
# Exercise Fluidsystem
-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
+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-fluidsystem a) as well as its mixture with water (exercise-fluidsystem b).
## Problem set-up
@@ -160,13 +160,18 @@ We now want to implement a pressure-dependent density for our component. Open th
$`\displaystyle \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 `2pproblem.hh` file by uncommenting line 54. Furthermore, the new component has to be set as a liquid phase in the fluid system, i.e. comment line 98 and uncomment line 99. The density distribution of this phase (rhoN) at the final simulation time should look like this:
+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 `2pproblem.hh` file by uncommenting line 54. Furthermore, the new component has to be set as a liquid phase in the fluid system, i.e. comment line 101 and uncomment line 102. The density distribution of this phase (rhoN) at the final simulation time should look like this:
+You can plot the density of the phase consisting of your compressible component by setting `PlotDensity` in `exercise-fluidsystem_a.input` to `true` and starting the simulation again.
+Compare the gnuplot output to the following plot of the density function from above:
+
+
+
### 3. Implement a new fluid system
-The problem file for this part of the exercise is `2p2cproblem.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_ (`ExerciseFluidsystemTwoPTwoCTypeTag`) from a _TypeTag_ (`TwoPTwoC`) that holds the miscible two-phase
+The problem file for this part of the exercise is `2p2cproblem.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_ (`ExerciseFluidsystemTwoPTwoCTypeTag`) from a _TypeTag_ (`TwoPTwoC`) that holds the miscible two-phase
two-component model properties:
```c++
diff --git a/exercises/exercise-fluidsystem/components/plotdensityfunction.py b/exercises/exercise-fluidsystem/components/plotdensityfunction.py
deleted file mode 100644
index 2b256689339b8c094d4bd3afb235be2b57fbb214..0000000000000000000000000000000000000000
--- a/exercises/exercise-fluidsystem/components/plotdensityfunction.py
+++ /dev/null
@@ -1,19 +0,0 @@
-#!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/exercises/exercise-fluidsystem/exercise-fluidsystem_a.input b/exercises/exercise-fluidsystem/exercise-fluidsystem_a.input
index c7c829043c00af4651a245b1b1da7036a6336ad1..8b18ca89524c9fefc63736d232498b67343b6939 100644
--- a/exercises/exercise-fluidsystem/exercise-fluidsystem_a.input
+++ b/exercises/exercise-fluidsystem/exercise-fluidsystem_a.input
@@ -8,3 +8,6 @@ Name = exercise-fluidsystem_a # name will be given to e.g. to the vtk result fil
[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
+
+[Output]
+PlotDensity = false # plot density over pressure for your component
diff --git a/exercises/extradoc/exercise-fluidsystem_a_densityfunction.png b/exercises/extradoc/exercise-fluidsystem_a_densityfunction.png
new file mode 100644
index 0000000000000000000000000000000000000000..8ea34fb8883516509d0215c3d5ff3dc379937c9d
Binary files /dev/null and b/exercises/extradoc/exercise-fluidsystem_a_densityfunction.png differ
diff --git a/exercises/solution/exercise-fluidsystem/2pproblem.hh b/exercises/solution/exercise-fluidsystem/2pproblem.hh
index 8a5ddee8e996c19192a68242f02f4530e9c71eaf..c228046e9de53e98f6611fc1d9346d1cf2c5c342 100644
--- a/exercises/solution/exercise-fluidsystem/2pproblem.hh
+++ b/exercises/solution/exercise-fluidsystem/2pproblem.hh
@@ -59,6 +59,9 @@
// The two-phase immiscible fluid system
#include <dumux/material/fluidsystems/2pimmiscible.hh>
+// The interface to create plots during simulation
+#include <dumux/io/gnuplotinterface.hh>
+
namespace Dumux{
// Forward declaration of the problem class
template <class TypeTag> class ExerciseFluidsystemProblemTwoP;
@@ -129,6 +132,7 @@ class ExerciseFluidsystemProblemTwoP : public PorousMediumFlowProblem<TypeTag>
using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry)::LocalView;
using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
+ using FluidState = typename GET_PROP_TYPE(TypeTag, FluidState);
enum {
waterPressureIdx = Indices::pressureIdx,
@@ -157,6 +161,10 @@ public:
// set the depth of the bottom of the reservoir
depthBOR_ = this->fvGridGeometry().bBoxMax()[dimWorld-1];
+
+ // plot density over pressure of the phase consisting of your component
+ if(getParam<bool>("Output.PlotDensity"))
+ plotDensity_();
}
/*!
@@ -283,8 +291,37 @@ public:
}
private:
+ void plotDensity_()
+ {
+ FluidState fluidState;
+ fluidState.setTemperature(temperature());
+ int numberOfPoints = 100;
+ Scalar xMin = 1e4;
+ Scalar xMax = 1e7;
+ Scalar spacing = std::pow((xMax/xMin), 1.0/(numberOfPoints-1));
+ std::vector<double> x(numberOfPoints);
+ std::vector<double> y(numberOfPoints);
+ for (int i=0; i<numberOfPoints; ++i)
+ x[i] = xMin*std::pow(spacing,i);
+ for (int i=0; i<numberOfPoints; ++i)
+ {
+ fluidState.setPressure(FluidSystem::phase1Idx, x[i]);
+ y[i] = FluidSystem::density(fluidState, FluidSystem::phase1Idx);
+ }
+ gnuplot_.resetPlot();
+ gnuplot_.setXRange(xMin, xMax);
+ gnuplot_.setOption("set logscale x 10");
+ gnuplot_.setOption("set key left top");
+ gnuplot_.setYRange(1440, 1480);
+ gnuplot_.setXlabel("pressure [Pa]");
+ gnuplot_.setYlabel("density [kg/m^3]");
+ gnuplot_.addDataSetToPlot(x, y, "YourComponentPhase_density.dat", "w l t 'Phase consisting of your component'");
+ gnuplot_.plot("YourComponentPhase_density");
+ }
+
Scalar eps_; //! small epsilon value
Scalar depthBOR_; //! depth at the bottom of the reservoir
+ Dumux::GnuplotInterface<double> gnuplot_; //! collects data for plotting
};
}
diff --git a/exercises/solution/exercise-fluidsystem/CMakeLists.txt b/exercises/solution/exercise-fluidsystem/CMakeLists.txt
index ff2ff0583be172aa0241158f59d912ecf56a6cc8..caccfb65baafc120fa7d628146be2de715749afa 100644
--- a/exercises/solution/exercise-fluidsystem/CMakeLists.txt
+++ b/exercises/solution/exercise-fluidsystem/CMakeLists.txt
@@ -3,13 +3,13 @@
dune_add_test(NAME exercise-fluidsystem_a_solution
SOURCES exercise-fluidsystem.cc
COMPILE_DEFINITIONS TYPETAG=ExerciseFluidsystemTwoPTypeTag
- CMD_ARGS exercise-fluidsystem_a.input)
+ CMD_ARGS exercise-fluidsystem_a_solution.input)
#part b: 2p2cproblem
dune_add_test(NAME exercise-fluidsystem_b_solution
SOURCES exercise-fluidsystem.cc
COMPILE_DEFINITIONS TYPETAG=ExerciseFluidsystemTwoPTwoCTypeTag
- CMD_ARGS exercise-fluidsystem_b.input)
+ CMD_ARGS exercise-fluidsystem_b_solution.input)
# add exercises to the common target
add_dependencies(test_exercises exercise-fluidsystem_a exercise-fluidsystem_b)
diff --git a/exercises/solution/exercise-fluidsystem/exercise-fluidsystem_a.input b/exercises/solution/exercise-fluidsystem/exercise-fluidsystem_a_solution.input
similarity index 81%
rename from exercises/solution/exercise-fluidsystem/exercise-fluidsystem_a.input
rename to exercises/solution/exercise-fluidsystem/exercise-fluidsystem_a_solution.input
index b137b0608e234f5c4f95bbac1e0c38a859602e9c..ccfa1908e5eff743991cc1b18771ebc6d58657d0 100644
--- a/exercises/solution/exercise-fluidsystem/exercise-fluidsystem_a.input
+++ b/exercises/solution/exercise-fluidsystem/exercise-fluidsystem_a_solution.input
@@ -8,3 +8,6 @@ Name = exercise-fluidsystem_a_solution # name will be given to e.g. to the vtk r
[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
+
+[Output]
+PlotDensity = true # plot density over pressure for your component
diff --git a/exercises/solution/exercise-fluidsystem/exercise-fluidsystem_b.input b/exercises/solution/exercise-fluidsystem/exercise-fluidsystem_b_solution.input
similarity index 100%
rename from exercises/solution/exercise-fluidsystem/exercise-fluidsystem_b.input
rename to exercises/solution/exercise-fluidsystem/exercise-fluidsystem_b_solution.input