Skip to content
Snippets Groups Projects
Commit 01dc15c4 authored by Beatrix Becker's avatar Beatrix Becker
Browse files

[exercise fluidsystem][2pproblem] add logarithmic spacing and get density from fluidsystem directly

parent d82bdbfe
No related branches found
No related tags found
1 merge request!21Feature/exercise fluidsystem gnuplot
......@@ -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
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment