From e395a5903743082a2de4612c44ffaf02aeec03f0 Mon Sep 17 00:00:00 2001
From: Beatrix Becker <beatrix.becker@iws.uni-stuttgart.de>
Date: Tue, 17 Jul 2018 14:25:20 +0200
Subject: [PATCH] [exercise fluidsystem][2pproblem] add gnuplot output also in
solution
---
.../exercise-fluidsystem/2pproblem.hh | 37 +++++++++++++++++++
.../exercise-fluidsystem_a.input | 3 ++
2 files changed, 40 insertions(+)
diff --git a/exercises/solution/exercise-fluidsystem/2pproblem.hh b/exercises/solution/exercise-fluidsystem/2pproblem.hh
index 8a5ddee8..c228046e 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/exercise-fluidsystem_a.input b/exercises/solution/exercise-fluidsystem/exercise-fluidsystem_a.input
index b137b060..ccfa1908 100644
--- a/exercises/solution/exercise-fluidsystem/exercise-fluidsystem_a.input
+++ b/exercises/solution/exercise-fluidsystem/exercise-fluidsystem_a.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
--
GitLab