diff --git a/exercises/exercise-fluidsystem/2pproblem.hh b/exercises/exercise-fluidsystem/2pproblem.hh index 0e5d7438aabeb190c60f5721457fd664b38a50a7..003047c007f3c02b7e82bf858a81e92f95210168 100644 --- a/exercises/exercise-fluidsystem/2pproblem.hh +++ b/exercises/exercise-fluidsystem/2pproblem.hh @@ -132,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, @@ -289,20 +290,30 @@ public: private: void plotDensity_() { - std::vector<double> x(100); - std::vector<double> y(100); - for (int i=0; i<100; ++i) - x[i] = 1e4 + i*100910; - for (int i=0; i<100; ++i) - y[i] = MyCompressibleComponent<Scalar>::liquidDensity(300, x[i]); + 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(1e4, 1e7); + 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, "YourComponent_density.dat", "w l t 'Liquid density of your component'"); - gnuplot_.plot("YourComponent_density"); + gnuplot_.addDataSetToPlot(x, y, "YourComponentPhase_density.dat", "w l t 'Phase consisting of your component'"); + gnuplot_.plot("YourComponentPhase_density"); } Scalar eps_; //! small epsilon value