Merge from upstream.
This commit is contained in:
commit
e502cc54d0
@ -50,6 +50,7 @@ opm/core/utility/cpgpreprocess/sparsetable.c \
|
||||
opm/core/utility/cpgpreprocess/facetopology.c \
|
||||
opm/core/utility/cpgpreprocess/uniquepoints.c \
|
||||
opm/core/utility/miscUtilities.cpp \
|
||||
opm/core/utility/SimulatorTimer.cpp \
|
||||
opm/core/utility/StopWatch.cpp \
|
||||
opm/core/utility/newwells.c \
|
||||
opm/core/utility/writeVtkData.cpp \
|
||||
@ -133,6 +134,7 @@ opm/core/utility/parameters/ParameterXML.hpp \
|
||||
opm/core/utility/parameters/tinyxml/tinystr.h \
|
||||
opm/core/utility/parameters/tinyxml/tinyxml.h \
|
||||
opm/core/utility/RootFinders.hpp \
|
||||
opm/core/utility/SimulatorTimer.hpp \
|
||||
opm/core/utility/SparseTable.hpp \
|
||||
opm/core/utility/SparseVector.hpp \
|
||||
opm/core/utility/StopWatch.hpp \
|
||||
|
@ -53,6 +53,7 @@
|
||||
#include <opm/core/WellsManager.hpp>
|
||||
#include <opm/core/newwells.h>
|
||||
#include <opm/core/utility/ErrorMacros.hpp>
|
||||
#include <opm/core/utility/SimulatorTimer.hpp>
|
||||
#include <opm/core/utility/StopWatch.hpp>
|
||||
#include <opm/core/utility/Units.hpp>
|
||||
#include <opm/core/utility/writeVtkData.hpp>
|
||||
@ -138,6 +139,37 @@ public:
|
||||
}
|
||||
}
|
||||
|
||||
// Initialize saturations so that there is water below woc,
|
||||
// and oil above.
|
||||
// TODO: add 'anitialiasing', obtaining a more precise woc
|
||||
// by f. ex. subdividing cells cut by the woc.
|
||||
void initWaterOilContact(const UnstructuredGrid& grid,
|
||||
const Opm::IncompPropertiesInterface& props,
|
||||
const double woc)
|
||||
{
|
||||
// Find out which cells should have water and which should have oil.
|
||||
std::vector<int> oil;
|
||||
std::vector<int> water;
|
||||
const int num_cells = grid.number_of_cells;
|
||||
oil.reserve(num_cells);
|
||||
water.reserve(num_cells);
|
||||
const int dim = grid.dimensions;
|
||||
for (int c = 0; c < num_cells; ++c) {
|
||||
const double z = grid.cell_centroids[dim*c + dim - 1];
|
||||
if (z > woc) {
|
||||
// Z is depth, we put water in the deepest parts
|
||||
// (even if oil is heavier...).
|
||||
water.push_back(c);
|
||||
} else {
|
||||
oil.push_back(c);
|
||||
}
|
||||
}
|
||||
|
||||
// Set saturations.
|
||||
setWaterSat(oil, props, MinSat);
|
||||
setWaterSat(water, props, MaxSat);
|
||||
}
|
||||
|
||||
int numPhases() const { return sat_.size()/press_.size(); }
|
||||
|
||||
::std::vector<double>& pressure () { return press_ ; }
|
||||
@ -159,11 +191,10 @@ private:
|
||||
|
||||
|
||||
|
||||
template <class State>
|
||||
void outputState(const UnstructuredGrid* grid,
|
||||
const State& state,
|
||||
const int step,
|
||||
const std::string& output_dir)
|
||||
void outputState(const UnstructuredGrid& grid,
|
||||
const ReservoirState& state,
|
||||
const int step,
|
||||
const std::string& output_dir)
|
||||
{
|
||||
// Write data in VTK format.
|
||||
std::ostringstream vtkfilename;
|
||||
@ -176,7 +207,7 @@ void outputState(const UnstructuredGrid* grid,
|
||||
dm["saturation"] = &state.saturation();
|
||||
dm["pressure"] = &state.pressure();
|
||||
std::vector<double> cell_velocity;
|
||||
Opm::estimateCellVelocity(*grid, state.faceflux(), cell_velocity);
|
||||
Opm::estimateCellVelocity(grid, state.faceflux(), cell_velocity);
|
||||
dm["velocity"] = &cell_velocity;
|
||||
Opm::writeVtkData(grid, dm, vtkfile);
|
||||
|
||||
@ -215,6 +246,39 @@ void wellsToSrc(const Wells& wells, const int num_cells, std::vector<double>& sr
|
||||
}
|
||||
}
|
||||
|
||||
class Watercut
|
||||
{
|
||||
public:
|
||||
void push(double time, double fraction, double produced)
|
||||
{
|
||||
data_.push_back(time);
|
||||
data_.push_back(fraction);
|
||||
data_.push_back(produced);
|
||||
}
|
||||
void write(std::ostream& os) const
|
||||
{
|
||||
int sz = data_.size()/3;
|
||||
for (int i = 0; i < sz; ++i) {
|
||||
os << data_[3*i]/Opm::unit::day << " "
|
||||
<< data_[3*i+1] << " "
|
||||
<< data_[3*i+2] << '\n';
|
||||
}
|
||||
}
|
||||
private:
|
||||
std::vector<double> data_;
|
||||
};
|
||||
|
||||
void outputWaterCut(const Watercut& watercut,
|
||||
const std::string& output_dir)
|
||||
{
|
||||
// Write water cut curve.
|
||||
std::string fname = output_dir + "/watercut.txt";
|
||||
std::ofstream os(fname.c_str());
|
||||
if (!os) {
|
||||
THROW("Failed to open " << fname);
|
||||
}
|
||||
watercut.write(os);
|
||||
}
|
||||
|
||||
// --------------- Types needed to define transport solver ---------------
|
||||
|
||||
@ -329,10 +393,7 @@ main(int argc, char** argv)
|
||||
Opm::parameter::ParameterGroup param(argc, argv, false);
|
||||
std::cout << "--------------- Reading parameters ---------------" << std::endl;
|
||||
|
||||
// Reading various control parameters.
|
||||
const int num_psteps = param.getDefault("num_psteps", 1);
|
||||
const double stepsize_days = param.getDefault("stepsize_days", 1.0);
|
||||
const double stepsize = Opm::unit::convert::from(stepsize_days, Opm::unit::day);
|
||||
// Reading various control parameters.
|
||||
const bool guess_old_solution = param.getDefault("guess_old_solution", false);
|
||||
const bool use_reorder = param.getDefault("use_reorder", true);
|
||||
const bool output = param.getDefault("output", true);
|
||||
@ -349,6 +410,9 @@ main(int argc, char** argv)
|
||||
boost::scoped_ptr<Opm::GridManager> grid;
|
||||
boost::scoped_ptr<Opm::IncompPropertiesInterface> props;
|
||||
boost::scoped_ptr<Opm::WellsManager> wells;
|
||||
Opm::SimulatorTimer simtimer;
|
||||
double water_oil_contact = 0.0;
|
||||
bool woc_set = false;
|
||||
if (use_deck) {
|
||||
std::string deck_filename = param.get<std::string>("deck_filename");
|
||||
Opm::EclipseGridParser deck(deck_filename);
|
||||
@ -360,6 +424,20 @@ main(int argc, char** argv)
|
||||
props.reset(new Opm::IncompPropertiesFromDeck(deck, global_cell));
|
||||
// Wells init.
|
||||
wells.reset(new Opm::WellsManager(deck, *grid->c_grid(), props->permeability()));
|
||||
// Timer init.
|
||||
if (deck.hasField("TSTEP")) {
|
||||
simtimer.init(deck);
|
||||
} else {
|
||||
simtimer.init(param);
|
||||
}
|
||||
// Water-oil contact.
|
||||
if (deck.hasField("EQUIL")) {
|
||||
water_oil_contact = deck.getEQUIL().equil[0].water_oil_contact_depth_;
|
||||
woc_set = true;
|
||||
} else if (param.has("water_oil_contact")) {
|
||||
water_oil_contact = param.get<double>("water_oil_contact");
|
||||
woc_set = true;
|
||||
}
|
||||
} else {
|
||||
// Grid init.
|
||||
const int nx = param.getDefault("nx", 100);
|
||||
@ -373,6 +451,12 @@ main(int argc, char** argv)
|
||||
props.reset(new Opm::IncompPropertiesBasic(param, grid->c_grid()->dimensions, grid->c_grid()->number_of_cells));
|
||||
// Wells init.
|
||||
wells.reset(new Opm::WellsManager());
|
||||
// Timer init.
|
||||
simtimer.init(param);
|
||||
if (param.has("water_oil_contact")) {
|
||||
water_oil_contact = param.get<double>("water_oil_contact");
|
||||
woc_set = true;
|
||||
}
|
||||
}
|
||||
|
||||
// Extra rock init.
|
||||
@ -445,11 +529,11 @@ main(int argc, char** argv)
|
||||
// code expects a scalar sw, not both sw and so.
|
||||
std::vector<double> reorder_sat(num_cells);
|
||||
std::vector<double> src(num_cells, 0.0);
|
||||
int scenario = param.getDefault("scenario", 0);
|
||||
int scenario = param.getDefault("scenario", woc_set ? 4 : 0);
|
||||
switch (scenario) {
|
||||
case 0:
|
||||
{
|
||||
std::cout << "==== Scenario 0: single-cell source and sink.\n";
|
||||
std::cout << "==== Scenario 0: simple wells or single-cell source and sink.\n";
|
||||
if (wells->c_wells()) {
|
||||
wellsToSrc(*wells->c_wells(), num_cells, src);
|
||||
} else {
|
||||
@ -477,6 +561,9 @@ main(int argc, char** argv)
|
||||
std::cout << "**** Warning: running gravity convection scenario, which expects a cartesian grid."
|
||||
<< std::endl;
|
||||
}
|
||||
if (grid->c_grid()->cartdims[2] <= 1) {
|
||||
std::cout << "**** Warning: running gravity convection scenario, which expects nz > 1." << std::endl;
|
||||
}
|
||||
std::vector<int> left_cells;
|
||||
left_cells.reserve(num_cells/2);
|
||||
const int *glob_cell = grid->c_grid()->global_cell;
|
||||
@ -501,6 +588,9 @@ main(int argc, char** argv)
|
||||
std::cout << "**** Warning: running gravity segregation scenario, which expects a cartesian grid."
|
||||
<< std::endl;
|
||||
}
|
||||
if (grid->c_grid()->cartdims[2] <= 1) {
|
||||
std::cout << "**** Warning: running gravity segregation scenario, which expects nz > 1." << std::endl;
|
||||
}
|
||||
std::vector<double>& sat = state.saturation();
|
||||
const int *glob_cell = grid->c_grid()->global_cell;
|
||||
// Water on top
|
||||
@ -513,6 +603,22 @@ main(int argc, char** argv)
|
||||
}
|
||||
break;
|
||||
}
|
||||
case 4:
|
||||
{
|
||||
std::cout << "==== Scenario 4: water-oil contact and simple wells or sources\n";
|
||||
if (!use_gravity) {
|
||||
std::cout << "**** Warning: initializing segregated water and oil zones, but gravity is zero." << std::endl;
|
||||
}
|
||||
state.initWaterOilContact(*grid->c_grid(), *props, water_oil_contact);
|
||||
if (wells->c_wells()) {
|
||||
wellsToSrc(*wells->c_wells(), num_cells, src);
|
||||
} else {
|
||||
double flow_per_sec = 0.01*tot_porevol/Opm::unit::day;
|
||||
src[0] = flow_per_sec;
|
||||
src[grid->c_grid()->number_of_cells - 1] = -flow_per_sec;
|
||||
}
|
||||
break;
|
||||
}
|
||||
default:
|
||||
{
|
||||
THROW("==== Scenario " << scenario << " is unknown.");
|
||||
@ -541,8 +647,6 @@ main(int argc, char** argv)
|
||||
// Control init.
|
||||
Opm::ImplicitTransportDetails::NRReport rpt;
|
||||
Opm::ImplicitTransportDetails::NRControl ctrl;
|
||||
double current_time = 0.0;
|
||||
double total_time = stepsize*num_psteps;
|
||||
if (!use_reorder || use_segregation_split) {
|
||||
ctrl.max_it = param.getDefault("max_it", 20);
|
||||
ctrl.verbosity = param.getDefault("verbosity", 0);
|
||||
@ -580,19 +684,20 @@ main(int argc, char** argv)
|
||||
Opm::time::StopWatch total_timer;
|
||||
total_timer.start();
|
||||
std::cout << "\n\n================ Starting main simulation loop ===============" << std::endl;
|
||||
double aver_sats[2] = { 0.0 };
|
||||
Opm::computeAverageSat(porevol, state.saturation(), aver_sats);
|
||||
std::cout << "\nInitial average saturations are " << aver_sats[0] << " " << aver_sats[1] << std::endl;
|
||||
for (int pstep = 0; pstep < num_psteps; ++pstep) {
|
||||
std::cout << "\n\n--------------- Simulation step number " << pstep
|
||||
<< " ---------------"
|
||||
<< "\n Current time (days) " << Opm::unit::convert::to(current_time, Opm::unit::day)
|
||||
<< "\n Current stepsize (days) " << Opm::unit::convert::to(stepsize, Opm::unit::day)
|
||||
<< "\n Total time (days) " << Opm::unit::convert::to(total_time, Opm::unit::day)
|
||||
<< "\n" << std::endl;
|
||||
|
||||
double satvol[2] = { 0.0 };
|
||||
double injected[2] = { 0.0 };
|
||||
double produced[2] = { 0.0 };
|
||||
double tot_injected[2] = { 0.0 };
|
||||
double tot_produced[2] = { 0.0 };
|
||||
Watercut watercut;
|
||||
watercut.push(0.0, 0.0, 0.0);
|
||||
Opm::computeSaturatedVol(porevol, state.saturation(), satvol);
|
||||
std::cout << "\nInitial saturated volumes are " << satvol[0] << " " << satvol[1] << std::endl;
|
||||
for (; !simtimer.done(); ++simtimer) {
|
||||
// Report timestep and (optionally) write state to disk.
|
||||
simtimer.report(std::cout);
|
||||
if (output) {
|
||||
outputState(grid->c_grid(), state, pstep, output_dir);
|
||||
outputState(*grid->c_grid(), state, simtimer.currentStepNum(), output_dir);
|
||||
}
|
||||
|
||||
// Solve pressure.
|
||||
@ -608,33 +713,39 @@ main(int argc, char** argv)
|
||||
std::cout << "Pressure solver took: " << pt << " seconds." << std::endl;
|
||||
ptime += pt;
|
||||
|
||||
// Process transport sources (to include bdy terms).
|
||||
if (use_reorder) {
|
||||
Opm::computeTransportSource(*grid->c_grid(), src, state.faceflux(), 1.0, reorder_src);
|
||||
} else {
|
||||
clear_transport_source(tsrc);
|
||||
for (int cell = 0; cell < num_cells; ++cell) {
|
||||
if (src[cell] > 0.0) {
|
||||
append_transport_source(cell, 2, 0, src[cell], ssrc, zdummy, tsrc);
|
||||
} else if (src[cell] < 0.0) {
|
||||
append_transport_source(cell, 2, 0, src[cell], ssink, zdummy, tsrc);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
// Solve transport
|
||||
transport_timer.start();
|
||||
if (use_reorder) {
|
||||
// We must treat reorder_src here,
|
||||
// if we are to handle anything but simple water
|
||||
// injection, since it is expected to be
|
||||
// equal to total outflow (if negative)
|
||||
// and water inflow (if positive).
|
||||
// Also, for anything but noflow boundaries,
|
||||
// boundary flows must be accumulated into
|
||||
// source term following the same convention.
|
||||
Opm::toWaterSat(state.saturation(), reorder_sat);
|
||||
reorder_model.solve(&state.faceflux()[0], &reorder_src[0], stepsize, &reorder_sat[0]);
|
||||
reorder_model.solve(&state.faceflux()[0], &reorder_src[0], simtimer.currentStepLength(), &reorder_sat[0]);
|
||||
Opm::toBothSat(reorder_sat, state.saturation());
|
||||
if (use_segregation_split) {
|
||||
if (use_column_solver) {
|
||||
colsolver.solve(columns, stepsize, state.saturation());
|
||||
colsolver.solve(columns, simtimer.currentStepLength(), state.saturation());
|
||||
} else {
|
||||
std::vector<double> fluxes = state.faceflux();
|
||||
std::fill(state.faceflux().begin(), state.faceflux().end(), 0.0);
|
||||
tsolver.solve(*grid->c_grid(), tsrc, stepsize, ctrl, state, linsolve, rpt);
|
||||
tsolver.solve(*grid->c_grid(), tsrc, simtimer.currentStepLength(), ctrl, state, linsolve, rpt);
|
||||
std::cout << rpt;
|
||||
state.faceflux() = fluxes;
|
||||
}
|
||||
}
|
||||
} else {
|
||||
tsolver.solve(*grid->c_grid(), tsrc, stepsize, ctrl, state, linsolve, rpt);
|
||||
tsolver.solve(*grid->c_grid(), tsrc, simtimer.currentStepLength(), ctrl, state, linsolve, rpt);
|
||||
std::cout << rpt;
|
||||
}
|
||||
transport_timer.stop();
|
||||
@ -642,10 +753,38 @@ main(int argc, char** argv)
|
||||
std::cout << "Transport solver took: " << tt << " seconds." << std::endl;
|
||||
ttime += tt;
|
||||
|
||||
Opm::computeAverageSat(porevol, state.saturation(), aver_sats);
|
||||
std::cout << "\nAverage saturations are " << aver_sats[0] << " " << aver_sats[1] << std::endl;
|
||||
Opm::computeSaturatedVol(porevol, state.saturation(), satvol);
|
||||
Opm::computeInjectedProduced(*props, state.saturation(), src, simtimer.currentStepLength(), injected, produced);
|
||||
tot_injected[0] += injected[0];
|
||||
tot_injected[1] += injected[1];
|
||||
tot_produced[0] += produced[0];
|
||||
tot_produced[1] += produced[1];
|
||||
std::cout.precision(5);
|
||||
const int width = 18;
|
||||
std::cout << "\nVolume balance report (all numbers relative to total pore volume).\n";
|
||||
std::cout << " Saturated volumes: "
|
||||
<< std::setw(width) << satvol[0]/tot_porevol
|
||||
<< std::setw(width) << satvol[1]/tot_porevol << std::endl;
|
||||
std::cout << " Injected volumes: "
|
||||
<< std::setw(width) << injected[0]/tot_porevol
|
||||
<< std::setw(width) << injected[1]/tot_porevol << std::endl;
|
||||
std::cout << " Produced volumes: "
|
||||
<< std::setw(width) << produced[0]/tot_porevol
|
||||
<< std::setw(width) << produced[1]/tot_porevol << std::endl;
|
||||
std::cout << " Total inj volumes: "
|
||||
<< std::setw(width) << tot_injected[0]/tot_porevol
|
||||
<< std::setw(width) << tot_injected[1]/tot_porevol << std::endl;
|
||||
std::cout << " Total prod volumes: "
|
||||
<< std::setw(width) << tot_produced[0]/tot_porevol
|
||||
<< std::setw(width) << tot_produced[1]/tot_porevol << std::endl;
|
||||
std::cout << " In-place + prod - inj: "
|
||||
<< std::setw(width) << (satvol[0] + tot_produced[0] - tot_injected[0])/tot_porevol
|
||||
<< std::setw(width) << (satvol[1] + tot_produced[1] - tot_injected[1])/tot_porevol << std::endl;
|
||||
std::cout.precision(8);
|
||||
|
||||
current_time += stepsize;
|
||||
watercut.push(simtimer.currentTime() + simtimer.currentStepLength(),
|
||||
produced[0]/(produced[0] + produced[1]),
|
||||
tot_produced[0]/tot_porevol);
|
||||
}
|
||||
total_timer.stop();
|
||||
|
||||
@ -655,7 +794,8 @@ main(int argc, char** argv)
|
||||
<< "\n Transport time: " << ttime << std::endl;
|
||||
|
||||
if (output) {
|
||||
outputState(grid->c_grid(), state, num_psteps, output_dir);
|
||||
outputState(*grid->c_grid(), state, simtimer.currentStepNum(), output_dir);
|
||||
outputWaterCut(watercut, output_dir);
|
||||
}
|
||||
|
||||
destroy_transport_source(tsrc);
|
||||
|
114
opm/core/utility/SimulatorTimer.cpp
Normal file
114
opm/core/utility/SimulatorTimer.cpp
Normal file
@ -0,0 +1,114 @@
|
||||
/*
|
||||
Copyright 2012 SINTEF ICT, Applied Mathematics.
|
||||
|
||||
This file is part of the Open Porous Media project (OPM).
|
||||
|
||||
OPM is free software: you can redistribute it and/or modify
|
||||
it under the terms of the GNU General Public License as published by
|
||||
the Free Software Foundation, either version 3 of the License, or
|
||||
(at your option) any later version.
|
||||
|
||||
OPM is distributed in the hope that it will be useful,
|
||||
but WITHOUT ANY WARRANTY; without even the implied warranty of
|
||||
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
||||
GNU General Public License for more details.
|
||||
|
||||
You should have received a copy of the GNU General Public License
|
||||
along with OPM. If not, see <http://www.gnu.org/licenses/>.
|
||||
*/
|
||||
|
||||
#include <opm/core/utility/SimulatorTimer.hpp>
|
||||
#include <opm/core/utility/parameters/ParameterGroup.hpp>
|
||||
#include <opm/core/utility/Units.hpp>
|
||||
#include <opm/core/eclipse/EclipseGridParser.hpp>
|
||||
#include <ostream>
|
||||
#include <numeric>
|
||||
|
||||
namespace Opm
|
||||
{
|
||||
|
||||
/// Default constructor.
|
||||
SimulatorTimer::SimulatorTimer()
|
||||
: current_step_(0),
|
||||
current_time_(0.0)
|
||||
{
|
||||
}
|
||||
|
||||
/// Initialize from parameters. Accepts the following:
|
||||
/// num_psteps (default 1)
|
||||
/// stepsize_days (default 1)
|
||||
void SimulatorTimer::init(const parameter::ParameterGroup& param)
|
||||
{
|
||||
const int num_psteps = param.getDefault("num_psteps", 1);
|
||||
const double stepsize_days = param.getDefault("stepsize_days", 1.0);
|
||||
const double stepsize = Opm::unit::convert::from(stepsize_days, Opm::unit::day);
|
||||
timesteps_.clear();
|
||||
timesteps_.resize(num_psteps, stepsize);
|
||||
total_time_ = num_psteps*stepsize;
|
||||
}
|
||||
|
||||
/// Initialize from TSTEP field.
|
||||
void SimulatorTimer::init(const EclipseGridParser& deck)
|
||||
{
|
||||
timesteps_ = deck.getTSTEP().tstep_;
|
||||
total_time_ = std::accumulate(timesteps_.begin(), timesteps_.end(), 0.0);
|
||||
}
|
||||
|
||||
/// Total number of steps.
|
||||
int SimulatorTimer::numSteps() const
|
||||
{
|
||||
return timesteps_.size();
|
||||
}
|
||||
|
||||
/// Current step number.
|
||||
int SimulatorTimer::currentStepNum() const
|
||||
{
|
||||
return current_step_;
|
||||
}
|
||||
|
||||
/// Current step length.
|
||||
double SimulatorTimer::currentStepLength() const
|
||||
{
|
||||
ASSERT(!done());
|
||||
return timesteps_[current_step_];
|
||||
}
|
||||
|
||||
/// Current time.
|
||||
double SimulatorTimer::currentTime() const
|
||||
{
|
||||
return current_time_;
|
||||
}
|
||||
|
||||
/// Total time.
|
||||
double SimulatorTimer::totalTime() const
|
||||
{
|
||||
return total_time_;
|
||||
}
|
||||
|
||||
/// Print a report with current and total time etc.
|
||||
void SimulatorTimer::report(std::ostream& os) const
|
||||
{
|
||||
os << "\n\n--------------- Simulation step number " << currentStepNum() << " ---------------"
|
||||
<< "\n Current time (days) " << Opm::unit::convert::to(currentTime(), Opm::unit::day)
|
||||
<< "\n Current stepsize (days) " << Opm::unit::convert::to(currentStepLength(), Opm::unit::day)
|
||||
<< "\n Total time (days) " << Opm::unit::convert::to(totalTime(), Opm::unit::day)
|
||||
<< "\n" << std::endl;
|
||||
}
|
||||
|
||||
/// Next step.
|
||||
SimulatorTimer& SimulatorTimer::operator++()
|
||||
{
|
||||
ASSERT(!done());
|
||||
current_time_ += timesteps_[current_step_];
|
||||
++current_step_;
|
||||
return *this;
|
||||
}
|
||||
|
||||
/// Return true if op++() has been called numSteps() times.
|
||||
bool SimulatorTimer::done() const
|
||||
{
|
||||
return int(timesteps_.size()) == current_step_;
|
||||
}
|
||||
|
||||
|
||||
} // namespace Opm
|
83
opm/core/utility/SimulatorTimer.hpp
Normal file
83
opm/core/utility/SimulatorTimer.hpp
Normal file
@ -0,0 +1,83 @@
|
||||
/*
|
||||
Copyright 2012 SINTEF ICT, Applied Mathematics.
|
||||
|
||||
This file is part of the Open Porous Media project (OPM).
|
||||
|
||||
OPM is free software: you can redistribute it and/or modify
|
||||
it under the terms of the GNU General Public License as published by
|
||||
the Free Software Foundation, either version 3 of the License, or
|
||||
(at your option) any later version.
|
||||
|
||||
OPM is distributed in the hope that it will be useful,
|
||||
but WITHOUT ANY WARRANTY; without even the implied warranty of
|
||||
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
||||
GNU General Public License for more details.
|
||||
|
||||
You should have received a copy of the GNU General Public License
|
||||
along with OPM. If not, see <http://www.gnu.org/licenses/>.
|
||||
*/
|
||||
|
||||
#ifndef OPM_SIMULATORTIMER_HEADER_INCLUDED
|
||||
#define OPM_SIMULATORTIMER_HEADER_INCLUDED
|
||||
|
||||
#include <iosfwd>
|
||||
#include <vector>
|
||||
|
||||
namespace Opm
|
||||
{
|
||||
|
||||
namespace parameter { class ParameterGroup; }
|
||||
class EclipseGridParser;
|
||||
|
||||
|
||||
class SimulatorTimer
|
||||
{
|
||||
public:
|
||||
/// Default constructor.
|
||||
SimulatorTimer();
|
||||
|
||||
/// Initialize from parameters. Accepts the following:
|
||||
/// num_psteps (default 1)
|
||||
/// stepsize_days (default 1)
|
||||
void init(const parameter::ParameterGroup& param);
|
||||
|
||||
/// Initialize from TSTEP field.
|
||||
void init(const EclipseGridParser& deck);
|
||||
|
||||
/// Total number of steps.
|
||||
int numSteps() const;
|
||||
|
||||
/// Current step number.
|
||||
int currentStepNum() const;
|
||||
|
||||
/// Current step length.
|
||||
/// Note: if done(), it is an error to call currentStepLength().
|
||||
double currentStepLength() const;
|
||||
|
||||
/// Current time.
|
||||
double currentTime() const;
|
||||
|
||||
/// Total time.
|
||||
double totalTime() const;
|
||||
|
||||
/// Print a report with current and total time etc.
|
||||
/// Note: if done(), it is an error to call report().
|
||||
void report(std::ostream& os) const;
|
||||
|
||||
/// Next step.
|
||||
SimulatorTimer& operator++();
|
||||
|
||||
/// Return true if op++() has been called numSteps() times.
|
||||
bool done() const;
|
||||
|
||||
private:
|
||||
std::vector<double> timesteps_;
|
||||
int current_step_;
|
||||
double current_time_;
|
||||
double total_time_;
|
||||
};
|
||||
|
||||
|
||||
} // namespace Opm
|
||||
|
||||
#endif // OPM_SIMULATORTIMER_HEADER_INCLUDED
|
@ -45,8 +45,33 @@ namespace Opm
|
||||
}
|
||||
|
||||
|
||||
/// @brief Computes total saturated volumes over all grid cells.
|
||||
/// @param[in] pv the pore volume by cell.
|
||||
/// @param[in] s saturation values (for all P phases)
|
||||
/// @param[out] sat_vol must point to a valid array with P elements,
|
||||
/// where P = s.size()/pv.size().
|
||||
/// For each phase p, we compute
|
||||
/// sat_vol_p = sum_i s_p_i pv_i
|
||||
void computeSaturatedVol(const std::vector<double>& pv,
|
||||
const std::vector<double>& s,
|
||||
double* sat_vol)
|
||||
{
|
||||
const int num_cells = pv.size();
|
||||
const int np = s.size()/pv.size();
|
||||
if (int(s.size()) != num_cells*np) {
|
||||
THROW("Sizes of s and pv vectors do not match.");
|
||||
}
|
||||
std::fill(sat_vol, sat_vol + np, 0.0);
|
||||
for (int c = 0; c < num_cells; ++c) {
|
||||
for (int p = 0; p < np; ++p) {
|
||||
sat_vol[p] += pv[c]*s[np*c + p];
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
/// @brief Computes average saturations over all grid cells.
|
||||
/// @param[out] pv the pore volume by cell.
|
||||
/// @param[in] pv the pore volume by cell.
|
||||
/// @param[in] s saturation values (for all P phases)
|
||||
/// @param[out] aver_sat must point to a valid array with P elements,
|
||||
/// where P = s.size()/pv.size().
|
||||
@ -78,6 +103,54 @@ namespace Opm
|
||||
}
|
||||
|
||||
|
||||
/// @brief Computes injected and produced volumes of all phases.
|
||||
/// Note 1: assumes that only the first phase is injected.
|
||||
/// Note 2: assumes that transport has been done with an
|
||||
/// implicit method, i.e. that the current state
|
||||
/// gives the mobilities used for the preceding timestep.
|
||||
/// @param[in] props fluid and rock properties.
|
||||
/// @param[in] s saturation values (for all P phases)
|
||||
/// @param[in] src if < 0: total outflow, if > 0: first phase inflow.
|
||||
/// @param[in] dt timestep used
|
||||
/// @param[out] injected must point to a valid array with P elements,
|
||||
/// where P = s.size()/src.size().
|
||||
/// @param[out] produced must also point to a valid array with P elements.
|
||||
void computeInjectedProduced(const IncompPropertiesInterface& props,
|
||||
const std::vector<double>& s,
|
||||
const std::vector<double>& src,
|
||||
const double dt,
|
||||
double* injected,
|
||||
double* produced)
|
||||
{
|
||||
const int num_cells = src.size();
|
||||
const int np = s.size()/src.size();
|
||||
if (int(s.size()) != num_cells*np) {
|
||||
THROW("Sizes of s and src vectors do not match.");
|
||||
}
|
||||
std::fill(injected, injected + np, 0.0);
|
||||
std::fill(produced, produced + np, 0.0);
|
||||
const double* visc = props.viscosity();
|
||||
std::vector<double> mob(np);
|
||||
for (int c = 0; c < num_cells; ++c) {
|
||||
if (src[c] > 0.0) {
|
||||
injected[0] += src[c]*dt;
|
||||
} else if (src[c] < 0.0) {
|
||||
const double flux = -src[c]*dt;
|
||||
const double* sat = &s[np*c];
|
||||
props.relperm(1, sat, &c, &mob[0], 0);
|
||||
double totmob = 0.0;
|
||||
for (int p = 0; p < np; ++p) {
|
||||
mob[p] /= visc[p];
|
||||
totmob += mob[p];
|
||||
}
|
||||
for (int p = 0; p < np; ++p) {
|
||||
produced[p] += (mob[p]/totmob)*flux;
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
|
||||
/// @brief Computes total mobility for a set of saturation values.
|
||||
/// @param[in] props rock and fluid properties
|
||||
@ -164,6 +237,48 @@ namespace Opm
|
||||
}
|
||||
}
|
||||
|
||||
/// Compute two-phase transport source terms from face fluxes,
|
||||
/// and pressure equation source terms. This puts boundary flows
|
||||
/// into the source terms for the transport equation.
|
||||
/// \param[in] grid The grid used.
|
||||
/// \param[in] src Pressure eq. source terms. The sign convention is:
|
||||
/// (+) positive total inflow (positive velocity divergence)
|
||||
/// (-) negative total outflow
|
||||
/// \param[in] faceflux Signed face fluxes, typically the result from a flow solver.
|
||||
/// \param[in] inflow_frac Fraction of inflow that consists of first phase.
|
||||
/// Example: if only water is injected, inflow_frac == 1.0.
|
||||
/// Note: it is not possible (with this method) to use different fractions
|
||||
/// for different inflow sources, be they source terms of boundary flows.
|
||||
/// \param[out] transport_src The transport source terms. They are to be interpreted depending on sign:
|
||||
/// (+) positive inflow of first phase (water)
|
||||
/// (-) negative total outflow of both phases
|
||||
void computeTransportSource(const UnstructuredGrid& grid,
|
||||
const std::vector<double>& src,
|
||||
const std::vector<double>& faceflux,
|
||||
const double inflow_frac,
|
||||
std::vector<double>& transport_src)
|
||||
{
|
||||
int nc = grid.number_of_cells;
|
||||
transport_src.resize(nc);
|
||||
for (int c = 0; c < nc; ++c) {
|
||||
transport_src[c] = 0.0;
|
||||
transport_src[c] += src[c] > 0.0 ? inflow_frac*src[c] : src[c];
|
||||
for (int hf = grid.cell_facepos[c]; hf < grid.cell_facepos[c + 1]; ++hf) {
|
||||
int f = grid.cell_faces[hf];
|
||||
const int* f2c = &grid.face_cells[2*f];
|
||||
double bdy_influx = 0.0;
|
||||
if (f2c[0] == c && f2c[1] == -1) {
|
||||
bdy_influx = -faceflux[f];
|
||||
} else if (f2c[0] == -1 && f2c[1] == c) {
|
||||
bdy_influx = faceflux[f];
|
||||
}
|
||||
if (bdy_influx != 0.0) {
|
||||
transport_src[c] += bdy_influx > 0.0 ? inflow_frac*bdy_influx : bdy_influx;
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
/// @brief Estimates a scalar cell velocity from face fluxes.
|
||||
/// @param[in] grid a grid
|
||||
/// @param[in] face_flux signed per-face fluxes
|
||||
|
@ -38,6 +38,18 @@ namespace Opm
|
||||
std::vector<double>& porevol);
|
||||
|
||||
|
||||
/// @brief Computes total saturated volumes over all grid cells.
|
||||
/// @param[out] pv the pore volume by cell.
|
||||
/// @param[in] s saturation values (for all P phases)
|
||||
/// @param[out] sat_vol must point to a valid array with P elements,
|
||||
/// where P = s.size()/pv.size().
|
||||
/// For each phase p, we compute
|
||||
/// sat_vol_p = sum_i s_p_i pv_i
|
||||
void computeSaturatedVol(const std::vector<double>& pv,
|
||||
const std::vector<double>& s,
|
||||
double* sat_vol);
|
||||
|
||||
|
||||
/// @brief Computes average saturations over all grid cells.
|
||||
/// @param[out] pv the pore volume by cell.
|
||||
/// @param[in] s saturation values (for all P phases)
|
||||
@ -50,6 +62,25 @@ namespace Opm
|
||||
double* aver_sat);
|
||||
|
||||
|
||||
/// @brief Computes injected and produced volumes of all phases.
|
||||
/// Note 1: assumes that only the first phase is injected.
|
||||
/// Note 2: assumes that transport has been done with an
|
||||
/// implicit method, i.e. that the current state
|
||||
/// gives the mobilities used for the preceding timestep.
|
||||
/// @param[in] props fluid and rock properties.
|
||||
/// @param[in] s saturation values (for all P phases)
|
||||
/// @param[in] src if < 0: total outflow, if > 0: first phase inflow.
|
||||
/// @param[in] dt timestep used
|
||||
/// @param[out] injected must point to a valid array with P elements,
|
||||
/// where P = s.size()/src.size().
|
||||
/// @param[out] produced must also point to a valid array with P elements.
|
||||
void computeInjectedProduced(const IncompPropertiesInterface& props,
|
||||
const std::vector<double>& s,
|
||||
const std::vector<double>& src,
|
||||
const double dt,
|
||||
double* injected,
|
||||
double* produced);
|
||||
|
||||
/// @brief Computes total mobility for a set of saturation values.
|
||||
/// @param[in] props rock and fluid properties
|
||||
/// @param[in] cells cells with which the saturation values are associated
|
||||
@ -73,11 +104,35 @@ namespace Opm
|
||||
std::vector<double>& totmob,
|
||||
std::vector<double>& omega);
|
||||
|
||||
|
||||
void computePhaseMobilities(const Opm::IncompPropertiesInterface& props,
|
||||
const std::vector<int>& cells,
|
||||
const std::vector<double>& s ,
|
||||
std::vector<double>& pmobc);
|
||||
|
||||
|
||||
/// Compute two-phase transport source terms from face fluxes,
|
||||
/// and pressure equation source terms. This puts boundary flows
|
||||
/// into the source terms for the transport equation.
|
||||
/// \param[in] grid The grid used.
|
||||
/// \param[in] src Pressure eq. source terms. The sign convention is:
|
||||
/// (+) positive total inflow (positive velocity divergence)
|
||||
/// (-) negative total outflow
|
||||
/// \param[in] faceflux Signed face fluxes, typically the result from a flow solver.
|
||||
/// \param[in] inflow_frac Fraction of inflow that consists of first phase.
|
||||
/// Example: if only water is injected, inflow_frac == 1.0.
|
||||
/// Note: it is not possible (with this method) to use different fractions
|
||||
/// for different inflow sources, be they source terms of boundary flows.
|
||||
/// \param[out] transport_src The transport source terms. They are to be interpreted depending on sign:
|
||||
/// (+) positive inflow of first phase (water)
|
||||
/// (-) negative total outflow of both phases
|
||||
void computeTransportSource(const UnstructuredGrid& grid,
|
||||
const std::vector<double>& src,
|
||||
const std::vector<double>& faceflux,
|
||||
const double inflow_frac,
|
||||
std::vector<double>& transport_src);
|
||||
|
||||
|
||||
/// @brief Estimates a scalar cell velocity from face fluxes.
|
||||
/// @param[in] grid a grid
|
||||
/// @param[in] face_flux signed per-face fluxes
|
||||
@ -91,7 +146,7 @@ namespace Opm
|
||||
void toWaterSat(const std::vector<double>& sboth,
|
||||
std::vector<double>& sw);
|
||||
|
||||
/// Make a a vector of interleaved water and oil saturations from
|
||||
/// Make a vector of interleaved water and oil saturations from
|
||||
/// a vector of water saturations.
|
||||
void toBothSat(const std::vector<double>& sw,
|
||||
std::vector<double>& sboth);
|
||||
|
@ -137,11 +137,11 @@ namespace Opm
|
||||
int Tag::indent_ = 0;
|
||||
|
||||
|
||||
void writeVtkData(const UnstructuredGrid* grid,
|
||||
void writeVtkData(const UnstructuredGrid& grid,
|
||||
const DataMap& data,
|
||||
std::ostream& os)
|
||||
{
|
||||
if (grid->dimensions != 3) {
|
||||
if (grid.dimensions != 3) {
|
||||
THROW("Vtk output for 3d grids only");
|
||||
}
|
||||
os.precision(12);
|
||||
@ -150,8 +150,8 @@ namespace Opm
|
||||
pm["type"] = "UnstructuredGrid";
|
||||
Tag vtkfiletag("VTKFile", pm, os);
|
||||
Tag ugtag("UnstructuredGrid", os);
|
||||
int num_pts = grid->number_of_nodes;
|
||||
int num_cells = grid->number_of_cells;
|
||||
int num_pts = grid.number_of_nodes;
|
||||
int num_cells = grid.number_of_cells;
|
||||
pm.clear();
|
||||
pm["NumberOfPoints"] = boost::lexical_cast<std::string>(num_pts);
|
||||
pm["NumberOfCells"] = boost::lexical_cast<std::string>(num_cells);
|
||||
@ -166,9 +166,9 @@ namespace Opm
|
||||
Tag datag("DataArray", pm, os);
|
||||
for (int i = 0; i < num_pts; ++i) {
|
||||
Tag::indent(os);
|
||||
os << grid->node_coordinates[3*i + 0] << ' '
|
||||
<< grid->node_coordinates[3*i + 1] << ' '
|
||||
<< grid->node_coordinates[3*i + 2] << '\n';
|
||||
os << grid.node_coordinates[3*i + 0] << ' '
|
||||
<< grid.node_coordinates[3*i + 1] << ' '
|
||||
<< grid.node_coordinates[3*i + 2] << '\n';
|
||||
}
|
||||
}
|
||||
{
|
||||
@ -185,10 +185,10 @@ namespace Opm
|
||||
int hf = 0;
|
||||
for (int c = 0; c < num_cells; ++c) {
|
||||
std::set<int> cell_pts;
|
||||
for (; hf < grid->cell_facepos[c+1]; ++hf) {
|
||||
int f = grid->cell_faces[hf];
|
||||
const int* fnbeg = grid->face_nodes + grid->face_nodepos[f];
|
||||
const int* fnend = grid->face_nodes + grid->face_nodepos[f+1];
|
||||
for (; hf < grid.cell_facepos[c+1]; ++hf) {
|
||||
int f = grid.cell_faces[hf];
|
||||
const int* fnbeg = grid.face_nodes + grid.face_nodepos[f];
|
||||
const int* fnend = grid.face_nodes + grid.face_nodepos[f+1];
|
||||
cell_pts.insert(fnbeg, fnend);
|
||||
}
|
||||
cell_numpts.push_back(cell_pts.size());
|
||||
@ -220,21 +220,21 @@ namespace Opm
|
||||
{
|
||||
pm["Name"] = "faces";
|
||||
Tag t("DataArray", pm, os);
|
||||
const int* fp = grid->cell_facepos;
|
||||
const int* fp = grid.cell_facepos;
|
||||
int offset = 0;
|
||||
for (int c = 0; c < num_cells; ++c) {
|
||||
Tag::indent(os);
|
||||
os << fp[c+1] - fp[c] << '\n';
|
||||
++offset;
|
||||
for (int hf = fp[c]; hf < fp[c+1]; ++hf) {
|
||||
int f = grid->cell_faces[hf];
|
||||
const int* np = grid->face_nodepos;
|
||||
int f = grid.cell_faces[hf];
|
||||
const int* np = grid.face_nodepos;
|
||||
int f_num_pts = np[f+1] - np[f];
|
||||
Tag::indent(os);
|
||||
os << f_num_pts << ' ';
|
||||
++offset;
|
||||
std::copy(grid->face_nodes + np[f],
|
||||
grid->face_nodes + np[f+1],
|
||||
std::copy(grid.face_nodes + np[f],
|
||||
grid.face_nodes + np[f+1],
|
||||
std::ostream_iterator<int>(os, " "));
|
||||
os << '\n';
|
||||
offset += f_num_pts;
|
||||
@ -290,7 +290,7 @@ namespace Opm
|
||||
for (DataMap::const_iterator dit = data.begin(); dit != data.end(); ++dit) {
|
||||
pm["Name"] = dit->first;
|
||||
const std::vector<double>& field = *(dit->second);
|
||||
const int num_comps = field.size()/grid->number_of_cells;
|
||||
const int num_comps = field.size()/grid.number_of_cells;
|
||||
pm["NumberOfComponents"] = boost::lexical_cast<std::string>(num_comps);
|
||||
Tag ptag("DataArray", pm, os);
|
||||
const int num_per_line = num_comps == 1 ? 5 : num_comps;
|
||||
|
@ -42,7 +42,7 @@ namespace Opm
|
||||
std::ostream& os);
|
||||
|
||||
/// Vtk output for general grids.
|
||||
void writeVtkData(const UnstructuredGrid* grid,
|
||||
void writeVtkData(const UnstructuredGrid& grid,
|
||||
const DataMap& data,
|
||||
std::ostream& os);
|
||||
} // namespace Opm
|
||||
|
Loading…
Reference in New Issue
Block a user