Commit 2a398553 authored by Beatrix Becker's avatar Beatrix Becker
Browse files

Merge branch 'feature/exercise-fluidsystem-gnuplot' into 'master'

Feature/exercise fluidsystem gnuplot

See merge request !21
parents d3655a5b 3d06620a
......@@ -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
......
......@@ -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
};
}
......
# 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:
![](../extradoc/exercise-fluidsystem_a_solution2.png)
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:
![](../extradoc/exercise-fluidsystem_a_densityfunction.png)
### 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++
......
#!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()
......@@ -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
......@@ -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
};
}
......
......@@ -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)
......
......@@ -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
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment