[exercise fluidsystem][2pproblem] add gnuplot output also in solution

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

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