Hi,
you want to plot the flux over a specific boundary over time, is that right?
You can do that in a postTimeStep-function, using either the
dumux/io/gnuplotinterface.hh or just writing the data in a file and
plotting with another program afterwards. In the postTimeStep you need
to iterate over all elements and all intersections until you find those
you want to calculate the flux over. Calculate and sum up the fluxes
over the interfaces and write them out for each time step.
I attach an example where this is done using the gnuplotinterface and
also an outputFile. You need to adapt it of course for your setup and
add the private variables:
Dumux::GnuplotInterface<double> gnuplot_;
std::ofstream outputFile_;
Also include the gnuplotinterface.hh.
Cheers
Beatrix
On 28.10.2015 05:36, Ait Mahiout Latifa wrote:
Hi,
is it possible to obtain curves for oil production at the output as a
function of time, after simulation with decoupled tutorial problem? If
yes, how we can obtain it?
Kind regards.
_______________________________________________
Dumux mailing list
[email protected]
https://listserv.uni-stuttgart.de/mailman/listinfo/dumux
void postTimeStep()
{
ParentType::postTimeStep();
//TODO get the velocity from the flux data (has to be calculated in the
pressure model)
Scalar inflow(0.0);
Scalar outflow(0.0);
Scalar simulationTime = this->timeManager().time();
Scalar temperature;
Scalar pressure;
// Calculate the inflow and outflow rate
for(auto eIt = this->gridView().template begin<0>(); eIt !=
this->gridView().template end<0>(); ++eIt)
{
int eIdx = this->variables().index(*eIt);
for(auto isIt = this->gridView().ibegin(*eIt); isIt !=
this->gridView().iend(*eIt); ++isIt)
{
if(isIt->boundary())
{
const GlobalPosition &globalPos =
isIt->geometry().center();
if(globalPos[2] < eps_)
{
outflow +=
(this->variables().cellData(eIdx).fluxData().velocity(wPhaseIdx,
isIt->indexInInside())
*isIt->centerUnitOuterNormal())*this->variables().cellData(eIdx).density(wPhaseIdx)
*isIt->geometry().volume();
if(globalPos[0] < eps_ && globalPos[1] < eps_)
temperature =
this->variables().cellData(eIdx).fluidState().temperature();
}
if(globalPos[2] > this->bBoxMax()[2] - eps_)
{
BoundaryTypes bcTypes;
this->boundaryTypes(bcTypes, *isIt);
if(globalPos[0] < eps_ && globalPos[1] < eps_)
pressure =
this->variables().cellData(eIdx).pressure(wPhaseIdx);
if(bcTypes.isNeumann(Indices::contiWEqIdx))
{
PrimaryVariables flux(0.0);
this->neumann(flux, *isIt);
inflow +=
flux[Indices::contiWEqIdx]*isIt->geometry().volume();
}
else
{
inflow +=
(this->variables().cellData(eIdx).fluxData().velocity(wPhaseIdx,
isIt->indexInInside())
*isIt->centerUnitOuterNormal())*this->variables().cellData(eIdx).density(wPhaseIdx)
*isIt->geometry().volume()*(-1.0);
}
}
}
}
}
std::cout << "Inflow rate: " << inflow*1000 << " g/s at t = " <<
simulationTime << std::endl;
std::cout << "Outflow rate: " << outflow*1000 << " g/s at t = " <<
simulationTime << std::endl;
outputFile_.open("flowdata.out", std::ios::app);
outputFile_ << simulationTime << " " << inflow*1000 << " " <<
outflow*1000 << "\n";
outputFile_.close();
// Gnuplot
// variables that don't go out of scope
static double yMax = -1e100;
static double y2Max = -1e100;
static std::vector<double> x;
static std::vector<double> y;
static std::vector<double> y2;
static std::vector<double> y3;
static std::vector<double> y4;
static std::vector<double> y5;
static double accumulatedMass = 0.0;
double inflow_plot = -1.0*inflow*1000;
double outflow_plot = -1.0*outflow*1000;
//PrimaryVariables bcValues;
//GlobalPosition globalPos({0.0, 0.0, this->bBoxMax()[2]});
//dirichletAtPos(bcValues, globalPos);
//double pressure = bcValues[Indices::pressureEqIdx]*1e-5;
pressure = pressure*1e-5;
temperature = temperature - 273.15;
accumulatedMass +=
outflow_plot*this->timeManager().previousTimeStepSize();
x.push_back(this->timeManager().time());
y.push_back(inflow_plot);
y2.push_back(pressure);
y3.push_back(accumulatedMass);
y4.push_back(temperature);
// y5.push_back(outflow_plot);
yMax = std::max({yMax, inflow_plot, pressure, outflow_plot});
y2Max = std::max({y2Max, accumulatedMass, temperature});
gnuplot_.reset();
gnuplot_.setXRange(0, x.back()+1);
gnuplot_.setYRange(0, yMax+1);
gnuplot_.setXlabel("time [s]");
gnuplot_.setYlabel("flowrate [g/s], pressure [bar]");
gnuplot_.addDataSetToPlot(x, y, "inflow at top", "axes x1y1 w l");
gnuplot_.addDataSetToPlot(x, y2, "pressure at top", "axes x1y1 w l");
gnuplot_.addDataSetToPlot(x, y3, "cumulated mass", "axes x1y2 w l");
gnuplot_.addDataSetToPlot(x, y4, "temperature at bottom", "axes x1y2 w
l");
// gnuplot_.addDataSetToPlot(x, y5, "outflow at bottom", "axes x1y1 w l");
// second axis
gnuplot_.setOption("set y2label \"temperature [°C], mass [g]\"");
std::ostringstream stream;
stream << "set y2range [" << 0 << ":" << y2Max+1 << "]";
gnuplot_.setOption(stream.str());
gnuplot_.setOption("set ytics nomirror");
gnuplot_.setOption("set y2tics");
gnuplot_.plot("inflow_plot");
}_______________________________________________
Dumux mailing list
[email protected]
https://listserv.uni-stuttgart.de/mailman/listinfo/dumux
