OilfieldToolsNet/opm-simulators
4
0
Fork 0
mirror of https://github.com/OPM/opm-simulators.git synced 2025-02-25 18:55:30 -06:00
Code Issues Packages Projects Releases Wiki Activity

merge

Browse Source
This commit is contained in:
Xavier Raynaud 2012-03-15 18:10:10 +01:00
commit dc0771ee3e
3 changed files with 661 additions and 249 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

748
examples/polymer_reorder.cpp
View File

@ -18,26 +18,33 @@
along with OPM. If not, see <http://www.gnu.org/licenses/>. along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/ */
#if HAVE_CONFIG_H
#include "config.h" #include "config.h"
#endif // HAVE_CONFIG_H
#include <opm/core/pressure/IncompTpfa.hpp> #include <opm/core/pressure/IncompTpfa.hpp>
#include <opm/core/pressure/FlowBCManager.hpp> #include <opm/core/pressure/FlowBCManager.hpp>
#include <opm/core/grid.h> #include <opm/core/grid.h>
#include <opm/core/GridManager.hpp> #include <opm/core/GridManager.hpp>
#include <opm/core/utility/writeVtkData.hpp> #include <opm/core/newwells.h>
#include <opm/core/utility/linearInterpolation.hpp> #include <opm/core/WellsManager.hpp>
#include <opm/core/utility/miscUtilities.hpp>
#include <opm/core/utility/ErrorMacros.hpp> #include <opm/core/utility/ErrorMacros.hpp>
#include <opm/core/utility/SimulatorTimer.hpp>
#include <opm/core/utility/StopWatch.hpp> #include <opm/core/utility/StopWatch.hpp>
#include <opm/core/utility/Units.hpp> #include <opm/core/utility/Units.hpp>
#include <opm/core/utility/writeVtkData.hpp>
#include <opm/core/utility/miscUtilities.hpp>
#include <opm/core/utility/parameters/ParameterGroup.hpp> #include <opm/core/utility/parameters/ParameterGroup.hpp>
#include <opm/core/fluid/IncompPropertiesBasic.hpp> #include <opm/core/fluid/IncompPropertiesBasic.hpp>
#include <opm/core/fluid/IncompPropertiesFromDeck.hpp> #include <opm/core/fluid/IncompPropertiesFromDeck.hpp>
#include <opm/core/linalg/LinearSolverUmfpack.hpp> #include <opm/core/linalg/LinearSolverUmfpack.hpp>
// #include <opm/core/linalg/LinearSolverIstl.hpp> // #define EXPERIMENT_ISTL
#ifdef EXPERIMENT_ISTL
#include <opm/core/linalg/LinearSolverIstl.hpp>
#endif
#include <opm/polymer/TransportModelPolymer.hpp> #include <opm/polymer/TransportModelPolymer.hpp>
#include <opm/polymer/PolymerProperties.hpp> #include <opm/polymer/PolymerProperties.hpp>
@ -66,79 +73,79 @@ class AdHocProps : public Opm::IncompPropertiesBasic
{ {
public: public:
AdHocProps(const Opm::parameter::ParameterGroup& param, int dim, int num_cells) AdHocProps(const Opm::parameter::ParameterGroup& param, int dim, int num_cells)
: Opm::IncompPropertiesBasic(param, dim, num_cells) : Opm::IncompPropertiesBasic(param, dim, num_cells)
{ {
ASSERT(numPhases() == 2); ASSERT(numPhases() == 2);
sw_.resize(3); sw_.resize(3);
sw_[0] = 0.2; sw_[0] = 0.2;
sw_[1] = 0.7; sw_[1] = 0.7;
sw_[2] = 1.0; sw_[2] = 1.0;
krw_.resize(3); krw_.resize(3);
krw_[0] = 0.0; krw_[0] = 0.0;
krw_[1] = 0.7; krw_[1] = 0.7;
krw_[2] = 1.0; krw_[2] = 1.0;
so_.resize(2); so_.resize(2);
so_[0] = 0.3; so_[0] = 0.3;
so_[1] = 0.8; so_[1] = 0.8;
kro_.resize(2); kro_.resize(2);
kro_[0] = 0.0; kro_[0] = 0.0;
kro_[1] = 1.0; kro_[1] = 1.0;
} }
virtual void relperm(const int n, virtual void relperm(const int n,
const double* s, const double* s,
const int* /*cells*/, const int* /*cells*/,
double* kr, double* kr,
double* dkrds) const double* dkrds) const
{ {
// ASSERT(dkrds == 0); // ASSERT(dkrds == 0);
// We assume two phases flow // We assume two phases flow
for (int i = 0; i < n; ++i) { for (int i = 0; i < n; ++i) {
kr[2*i] = krw(s[2*i]); kr[2*i] = krw(s[2*i]);
kr[2*i+1] = kro(s[2*i+1]); kr[2*i+1] = kro(s[2*i+1]);
if (dkrds != 0) { if (dkrds != 0) {
dkrds[2*i] = krw_dsw(s[2*i]); dkrds[2*i] = krw_dsw(s[2*i]);
dkrds[2*i+3] = kro_dso(s[2*i+1]); dkrds[2*i+3] = kro_dso(s[2*i+1]);
dkrds[2*i+1] = -dkrds[2*i+3]; dkrds[2*i+1] = -dkrds[2*i+3];
dkrds[2*i+2] = -dkrds[2*i]; dkrds[2*i+2] = -dkrds[2*i];
} }
} }
} }
virtual void satRange(const int n, virtual void satRange(const int n,
const int* /*cells*/, const int* /*cells*/,
double* smin, double* smin,
double* smax) const double* smax) const
{ {
const int np = 2; const int np = 2;
for (int i = 0; i < n; ++i) { for (int i = 0; i < n; ++i) {
smin[np*i + 0] = sw_[0]; smin[np*i + 0] = sw_[0];
smax[np*i + 0] = sw_.back(); smax[np*i + 0] = sw_.back();
smin[np*i + 1] = 1.0 - sw_[0]; smin[np*i + 1] = 1.0 - sw_[0];
smax[np*i + 1] = 1.0 - sw_.back(); smax[np*i + 1] = 1.0 - sw_.back();
} }
} }
private: private:
double krw(double s) const double krw(double s) const
{ {
return Opm::linearInterpolation(sw_, krw_, s); return Opm::linearInterpolation(sw_, krw_, s);
} }
double krw_dsw(double s) const double krw_dsw(double s) const
{ {
return Opm::linearInterpolationDerivative(sw_, krw_, s); return Opm::linearInterpolationDerivative(sw_, krw_, s);
} }
double kro(double s) const double kro(double s) const
{ {
return Opm::linearInterpolation(so_, kro_, s); return Opm::linearInterpolation(so_, kro_, s);
} }
double kro_dso(double s) const double kro_dso(double s) const
{ {
return Opm::linearInterpolationDerivative(so_, kro_, s); return Opm::linearInterpolationDerivative(so_, kro_, s);
} }
std::vector<double> sw_; std::vector<double> sw_;
@ -156,29 +163,72 @@ public:
fpress_(g->number_of_faces, 0.0), fpress_(g->number_of_faces, 0.0),
flux_ (g->number_of_faces, 0.0), flux_ (g->number_of_faces, 0.0),
sat_ (2 * g->number_of_cells, 0.0), sat_ (2 * g->number_of_cells, 0.0),
concentration_(g->number_of_cells, 0.0), concentration_(g->number_of_cells, 0.0),
cmax_(g->number_of_cells, 0.0) cmax_(g->number_of_cells, 0.0)
{ {
for (int cell = 0; cell < g->number_of_cells; ++cell) { for (int cell = 0; cell < g->number_of_cells; ++cell) {
sat_[2*cell] = init_sat; sat_[2*cell] = init_sat;
sat_[2*cell + 1] = 1.0 - init_sat; sat_[2*cell + 1] = 1.0 - init_sat;
} }
} }
enum ExtremalSat { MinSat, MaxSat };
void setToMinimumWaterSat(const Opm::IncompPropertiesInterface& props) void setToMinimumWaterSat(const Opm::IncompPropertiesInterface& props)
{ {
const int n = props.numCells(); const int n = props.numCells();
std::vector<int> cells(n); std::vector<int> cells(n);
for (int i = 0; i < n; ++i) { for (int i = 0; i < n; ++i) {
cells[i] = i; cells[i] = i;
} }
std::vector<double> smin(2*n); setWaterSat(cells, props, MinSat);
std::vector<double> smax(2*n); }
props.satRange(n, &cells[0], &smin[0], &smax[0]);
for (int cell = 0; cell < n; ++cell) { void setWaterSat(const std::vector<int>& cells,
sat_[2*cell] = smin[2*cell]; const Opm::IncompPropertiesInterface& props,
sat_[2*cell + 1] = 1.0 - smin[2*cell]; ExtremalSat es)
} {
const int n = cells.size();
std::vector<double> smin(2*n);
std::vector<double> smax(2*n);
props.satRange(n, &cells[0], &smin[0], &smax[0]);
const double* svals = (es == MinSat) ? &smin[0] : &smax[0];
for (int ci = 0; ci < n; ++ci) {
const int cell = cells[ci];
sat_[2*cell] = svals[2*ci];
sat_[2*cell + 1] = 1.0 - sat_[2*cell];
}
}
// 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(); } int numPhases() const { return sat_.size()/press_.size(); }
@ -208,22 +258,61 @@ private:
static void outputState(const UnstructuredGrid& grid,
const ReservoirState& state,
const int step,
const std::string& output_dir)
{
// Write data in VTK format.
std::ostringstream vtkfilename;
vtkfilename << output_dir << "/output-" << std::setw(3) << std::setfill('0') << step << ".vtu";
std::ofstream vtkfile(vtkfilename.str().c_str());
if (!vtkfile) {
THROW("Failed to open " << vtkfilename.str());
}
Opm::DataMap dm;
dm["saturation"] = &state.saturation();
dm["pressure"] = &state.pressure();
dm["concentration"] = &state.concentration();
dm["cmax"] = &state.cmax();
std::vector<double> cell_velocity;
Opm::estimateCellVelocity(grid, state.faceflux(), cell_velocity);
dm["velocity"] = &cell_velocity;
Opm::writeVtkData(grid, dm, vtkfile);
// Write data (not grid) in Matlab format
for (Opm::DataMap::const_iterator it = dm.begin(); it != dm.end(); ++it) {
std::ostringstream fname;
fname << output_dir << "/" << it->first << "-" << std::setw(3) << std::setfill('0') << step << ".dat";
std::ofstream file(fname.str().c_str());
if (!file) {
THROW("Failed to open " << fname.str());
}
const std::vector<double>& d = *(it->second);
std::copy(d.begin(), d.end(), std::ostream_iterator<double>(file, "\n"));
}
}
class PolymerInflow class PolymerInflow
{ {
public: public:
PolymerInflow(const double starttime, PolymerInflow(const double starttime,
const double endtime, const double endtime,
const double amount) const double amount)
: stime_(starttime), etime_(endtime), amount_(amount) : stime_(starttime), etime_(endtime), amount_(amount)
{ {
} }
double operator()(double time) double operator()(double time)
{ {
if (time >= stime_ && time < etime_) { if (time >= stime_ && time < etime_) {
return amount_; return amount_;
} else { } else {
return 0.0; return 0.0;
} }
} }
private: private:
double stime_; double stime_;
@ -233,44 +322,19 @@ private:
static void outputWaterCut(const Opm::Watercut& watercut,
template <class State> const std::string& output_dir)
void outputState(const UnstructuredGrid* grid,
const State& state,
const int step,
const std::string& output_dir)
{ {
// Write data in VTK format. // Write water cut curve.
std::ostringstream vtkfilename; std::string fname = output_dir + "/watercut.txt";
vtkfilename << output_dir << "/output-" << std::setw(3) << std::setfill('0') << step << ".vtu"; std::ofstream os(fname.c_str());
std::ofstream vtkfile(vtkfilename.str().c_str()); if (!os) {
if (!vtkfile) { THROW("Failed to open " << fname);
THROW("Failed to open " << vtkfilename.str());
}
Opm::DataMap dm;
dm["saturation"] = &state.saturation();
dm["pressure"] = &state.pressure();
dm["concentration"] = &state.concentration();
Opm::writeVtkData(grid, dm, vtkfile);
// Write data (not grid) in Matlab format
for (Opm::DataMap::const_iterator it = dm.begin(); it != dm.end(); ++it) {
std::ostringstream fname;
fname << output_dir << "/" << it->first << "-" << std::setw(3) << std::setfill('0') << step << ".dat";
std::ofstream file(fname.str().c_str());
if (!file) {
THROW("Failed to open " << fname.str());
}
const std::vector<double>& d = *(it->second);
std::copy(d.begin(), d.end(), std::ostream_iterator<double>(file, "\n"));
} }
watercut.write(os);
} }
// ----------------- Main program ----------------- // ----------------- Main program -----------------
int int
main(int argc, char** argv) main(int argc, char** argv)
@ -280,83 +344,101 @@ main(int argc, char** argv)
std::cout << "--------------- Reading parameters ---------------" << std::endl; std::cout << "--------------- Reading parameters ---------------" << std::endl;
// Reading various control parameters. // 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);
const bool output = param.getDefault("output", true); const bool output = param.getDefault("output", true);
std::string output_dir; std::string output_dir;
int output_interval = 1;
if (output) { if (output) {
output_dir = param.getDefault("output_dir", std::string("output")); output_dir = param.getDefault("output_dir", std::string("output"));
// Ensure that output dir exists // Ensure that output dir exists
boost::filesystem::path fpath(output_dir); boost::filesystem::path fpath(output_dir);
create_directories(fpath); create_directories(fpath);
output_interval = param.getDefault("output_interval", output_interval);
} }
// If we have a "deck_filename", grid and props will be read from that. // If we have a "deck_filename", grid and props will be read from that.
bool use_deck = param.has("deck_filename"); bool use_deck = param.has("deck_filename");
boost::scoped_ptr<Opm::GridManager> grid; boost::scoped_ptr<Opm::GridManager> grid;
boost::scoped_ptr<Opm::IncompPropertiesInterface> props; 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;
Opm::PolymerProperties polydata; Opm::PolymerProperties polydata;
if (use_deck) { if (use_deck) {
std::string deck_filename = param.get<std::string>("deck_filename"); std::string deck_filename = param.get<std::string>("deck_filename");
Opm::EclipseGridParser deck(deck_filename); Opm::EclipseGridParser deck(deck_filename);
polydata.readFromDeck(deck); // Grid init
// Grid init grid.reset(new Opm::GridManager(deck));
// grid.reset(new Opm::GridManager(deck)); // Rock and fluid init
const int nx = param.getDefault("nx", 100); const int* gc = grid->c_grid()->global_cell;
const int ny = param.getDefault("ny", 100); std::vector<int> global_cell(gc, gc + grid->c_grid()->number_of_cells);
const int nz = param.getDefault("nz", 1); props.reset(new Opm::IncompPropertiesFromDeck(deck, global_cell));
const double dx = param.getDefault("dx", 1.0); // props.reset(new AdHocProps(param, grid->c_grid()->dimensions, grid->c_grid()->number_of_cells));
const double dy = param.getDefault("dy", 1.0); // Wells init.
const double dz = param.getDefault("dz", 1.0); wells.reset(new Opm::WellsManager(deck, *grid->c_grid(), props->permeability()));
grid.reset(new Opm::GridManager(nx, ny, nz, dx, dy, dz)); // Timer init.
// Rock and fluid init if (deck.hasField("TSTEP")) {
const int* gc = grid->c_grid()->global_cell; simtimer.init(deck);
std::vector<int> global_cell(gc, gc + grid->c_grid()->number_of_cells); } else {
props.reset(new Opm::IncompPropertiesFromDeck(deck, global_cell)); simtimer.init(param);
// props.reset(new AdHocProps(param, grid->c_grid()->dimensions, grid->c_grid()->number_of_cells)); }
// 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;
}
polydata.readFromDeck(deck);
} else { } else {
// Grid init. // Grid init.
const int nx = param.getDefault("nx", 100); const int nx = param.getDefault("nx", 100);
const int ny = param.getDefault("ny", 100); const int ny = param.getDefault("ny", 100);
const int nz = param.getDefault("nz", 1); const int nz = param.getDefault("nz", 1);
const double dx = param.getDefault("dx", 1.0); const double dx = param.getDefault("dx", 1.0);
const double dy = param.getDefault("dy", 1.0); const double dy = param.getDefault("dy", 1.0);
const double dz = param.getDefault("dz", 1.0); const double dz = param.getDefault("dz", 1.0);
grid.reset(new Opm::GridManager(nx, ny, nz, dx, dy, dz)); grid.reset(new Opm::GridManager(nx, ny, nz, dx, dy, dz));
// Rock and fluid init. // Rock and fluid init.
// props.reset(new Opm::IncompPropertiesBasic(param, grid->c_grid()->dimensions, grid->c_grid()->number_of_cells)); // props.reset(new Opm::IncompPropertiesBasic(param, grid->c_grid()->dimensions, grid->c_grid()->number_of_cells));
props.reset(new AdHocProps(param, grid->c_grid()->dimensions, grid->c_grid()->number_of_cells)); props.reset(new AdHocProps(param, grid->c_grid()->dimensions, grid->c_grid()->number_of_cells));
// Setting polydata defaults to mimic a simple example case. // Wells init.
wells.reset(new Opm::WellsManager());
double c_max = param.getDefault("c_max_limit", 5.0); // Timer init.
double mix_param = param.getDefault("mix_param", 1.0); simtimer.init(param);
double rock_density = param.getDefault("rock_density", 1000.0); if (param.has("water_oil_contact")) {
double dead_pore_vol = param.getDefault("dead_pore_vol", 0.15); water_oil_contact = param.get<double>("water_oil_contact");
std::vector<double> c_vals_visc(2, -1e100); woc_set = true;
c_vals_visc[0] = 0.0; }
c_vals_visc[1] = 7.0; // Setting polydata defaults to mimic a simple example case.
std::vector<double> visc_mult_vals(2, -1e100); double c_max = param.getDefault("c_max_limit", 5.0);
visc_mult_vals[0] = 1.0; double mix_param = param.getDefault("mix_param", 1.0);
// polydata.visc_mult_vals[1] = param.getDefault("c_max_viscmult", 30.0); double rock_density = param.getDefault("rock_density", 1000.0);
visc_mult_vals[1] = 20.0; double dead_pore_vol = param.getDefault("dead_pore_vol", 0.15);
std::vector<double> c_vals_ads(3, -1e100); std::vector<double> c_vals_visc(2, -1e100);
c_vals_ads[0] = 0.0; c_vals_visc[0] = 0.0;
c_vals_ads[1] = 2.0; c_vals_visc[1] = 7.0;
c_vals_ads[2] = 8.0; std::vector<double> visc_mult_vals(2, -1e100);
std::vector<double> ads_vals(3, -1e100); visc_mult_vals[0] = 1.0;
ads_vals[0] = 0.0; // polydata.visc_mult_vals[1] = param.getDefault("c_max_viscmult", 30.0);
// polydata.ads_vals[1] = param.getDefault("c_max_ads", 0.0025); visc_mult_vals[1] = 20.0;
ads_vals[1] = 0.0015; std::vector<double> c_vals_ads(3, -1e100);
ads_vals[2] = 0.0025; c_vals_ads[0] = 0.0;
polydata.set(c_max, mix_param, rock_density, dead_pore_vol, c_vals_visc, visc_mult_vals, c_vals_ads, ads_vals); c_vals_ads[1] = 2.0;
c_vals_ads[2] = 8.0;
std::vector<double> ads_vals(3, -1e100);
ads_vals[0] = 0.0;
// polydata.ads_vals[1] = param.getDefault("c_max_ads", 0.0025);
ads_vals[1] = 0.0015;
ads_vals[2] = 0.0025;
polydata.set(c_max, mix_param, rock_density, dead_pore_vol, c_vals_visc, visc_mult_vals, c_vals_ads, ads_vals);
} }
// Initialize polymer inflow function.
double poly_start = param.getDefault("poly_start_days", 300.0)*Opm::unit::day; double poly_start = param.getDefault("poly_start_days", 300.0)*Opm::unit::day;
double poly_end = param.getDefault("poly_end_days", 800.0)*Opm::unit::day; double poly_end = param.getDefault("poly_end_days", 800.0)*Opm::unit::day;
double poly_amount = param.getDefault("poly_amount", 5.0); double poly_amount = param.getDefault("poly_amount", polydata.cMax());
PolymerInflow poly_inflow(poly_start, poly_end, poly_amount); PolymerInflow poly_inflow(poly_start, poly_end, poly_amount);
// Extra rock init. // Extra rock init.
@ -369,31 +451,35 @@ main(int argc, char** argv)
double g = param.getDefault("gravity", 0.0); double g = param.getDefault("gravity", 0.0);
bool use_gravity = g != 0.0; bool use_gravity = g != 0.0;
if (use_gravity) { if (use_gravity) {
gravity[grid->c_grid()->dimensions - 1] = g; gravity[grid->c_grid()->dimensions - 1] = g;
if (props->density()[0] == props->density()[1]) { if (props->density()[0] == props->density()[1]) {
std::cout << "**** Warning: nonzero gravity, but zero density difference." << std::endl; std::cout << "**** Warning: nonzero gravity, but zero density difference." << std::endl;
} }
} }
// Solvers init. // Solvers init.
// Pressure solver.
#ifdef EXPERIMENT_ISTL
Opm::LinearSolverIstl linsolver(param);
#else
Opm::LinearSolverUmfpack linsolver; Opm::LinearSolverUmfpack linsolver;
// Opm::LinearSolverIstl linsolver(param); #endif // EXPERIMENT_ISTL
const double *grav = use_gravity ? &gravity[0] : 0; const double *grav = use_gravity ? &gravity[0] : 0;
Opm::IncompTpfa psolver(*grid->c_grid(), props->permeability(), grav, linsolver); Opm::IncompTpfa psolver(*grid->c_grid(), props->permeability(), grav, linsolver);
// Reordering solver.
const double nltol = param.getDefault("nl_tolerance", 1e-9);
const int maxit = param.getDefault("nl_maxiter", 30);
Opm::TransportModelPolymer::SingleCellMethod method; Opm::TransportModelPolymer::SingleCellMethod method;
std::string method_string = param.getDefault("single_cell_method", std::string("Bracketing")); std::string method_string = param.getDefault("single_cell_method", std::string("Bracketing"));
if (method_string == "Bracketing") { if (method_string == "Bracketing") {
method = Opm::TransportModelPolymer::Bracketing; method = Opm::TransportModelPolymer::Bracketing;
} else if (method_string == "Newton") { } else if (method_string == "Newton") {
method = Opm::TransportModelPolymer::Newton; method = Opm::TransportModelPolymer::Newton;
} else { } else {
THROW("Unknown method: " << method_string); THROW("Unknown method: " << method_string);
} }
const double nltol = param.getDefault("nl_tolerance", 1e-9);
const int maxit = param.getDefault("nl_maxiter", 30);
Opm::TransportModelPolymer tmodel(*grid->c_grid(), props->porosity(), &porevol[0], *props, polydata, Opm::TransportModelPolymer tmodel(*grid->c_grid(), props->porosity(), &porevol[0], *props, polydata,
method, nltol, maxit); method, nltol, maxit);
// Boundary conditions. // Boundary conditions.
Opm::FlowBCManager bcs; Opm::FlowBCManager bcs;
@ -401,45 +487,140 @@ main(int argc, char** argv)
// State-related and source-related variables init. // State-related and source-related variables init.
int num_cells = grid->c_grid()->number_of_cells; int num_cells = grid->c_grid()->number_of_cells;
std::vector<double> totmob; std::vector<double> totmob;
std::vector<double> omega; // Empty dummy unless/until we include gravity here. std::vector<double> omega; // Will remain empty if no gravity.
double init_sat = param.getDefault("init_sat", 0.0); double init_sat = param.getDefault("init_sat", 0.0);
ReservoirState state(grid->c_grid(), init_sat); ReservoirState state(grid->c_grid(), init_sat);
if (!param.has("init_sat")) { if (!param.has("init_sat")) {
state.setToMinimumWaterSat(*props); state.setToMinimumWaterSat(*props);
} }
// We need a separate reorder_sat, because the reorder // We need a separate reorder_sat, because the reorder
// code expects a scalar sw, not both sw and so. // code expects a scalar sw, not both sw and so.
std::vector<double> reorder_sat(num_cells); std::vector<double> reorder_sat(num_cells);
double flow_per_sec = 0.1*tot_porevol/Opm::unit::day;
if (param.has("injection_rate_per_day")) {
flow_per_sec = param.get<double>("injection_rate_per_day")/Opm::unit::day;
}
std::vector<double> src(num_cells, 0.0); std::vector<double> src(num_cells, 0.0);
src[0] = flow_per_sec; int scenario = param.getDefault("scenario", woc_set ? 4 : 0);
src[num_cells - 1] = -flow_per_sec; switch (scenario) {
case 0:
{
std::cout << "==== Scenario 0: simple wells or single-cell source and sink.\n";
if (wells->c_wells()) {
Opm::wellsToSrc(*wells->c_wells(), num_cells, src);
} else {
double flow_per_sec = 0.1*tot_porevol/Opm::unit::day;
if (param.has("injection_rate_per_day")) {
flow_per_sec = param.get<double>("injection_rate_per_day")/Opm::unit::day;
}
src[0] = flow_per_sec;
src[num_cells - 1] = -flow_per_sec;
}
break;
}
case 1:
{
std::cout << "==== Scenario 1: half source, half sink.\n";
double flow_per_sec = 0.1*porevol[0]/Opm::unit::day;
std::fill(src.begin(), src.begin() + src.size()/2, flow_per_sec);
std::fill(src.begin() + src.size()/2, src.end(), -flow_per_sec);
break;
}
case 2:
{
std::cout << "==== Scenario 2: gravity convection.\n";
if (!use_gravity) {
std::cout << "**** Warning: running gravity convection scenario, but gravity is zero." << std::endl;
}
if (use_deck) {
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;
for (int cell = 0; cell < num_cells; ++cell) {
const int* cd = grid->c_grid()->cartdims;
const int gc = glob_cell == 0 ? cell : glob_cell[cell];
bool left = (gc % cd[0]) < cd[0]/2;
if (left) {
left_cells.push_back(cell);
}
}
state.setWaterSat(left_cells, *props, ReservoirState::MaxSat);
break;
}
case 3:
{
std::cout << "==== Scenario 3: gravity segregation.\n";
if (!use_gravity) {
std::cout << "**** Warning: running gravity segregation scenario, but gravity is zero." << std::endl;
}
if (use_deck) {
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
for (int cell = 0; cell < num_cells; ++cell) {
const int* cd = grid->c_grid()->cartdims;
const int gc = glob_cell == 0 ? cell : glob_cell[cell];
bool top = (gc / cd[0] / cd[1]) < cd[2]/2;
sat[2*cell] = top ? 1.0 : 0.0;
sat[2*cell + 1 ] = 1.0 - sat[2*cell];
}
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()) {
Opm::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.");
}
}
std::vector<double> reorder_src = src; std::vector<double> reorder_src = src;
// Control init. // Dirichlet boundary conditions.
double current_time = 0.0; if (param.getDefault("use_pside", false)) {
double total_time = stepsize*num_psteps; int pside = param.get<int>("pside");
double pside_pressure = param.get<double>("pside_pressure");
bcs.pressureSide(*grid->c_grid(), Opm::FlowBCManager::Side(pside), pside_pressure);
}
// The allcells vector is used in calls to computeTotalMobility() // The allcells vector is used in calls to computeTotalMobility()
// and computeTotalMobilityOmega(). // and computeTotalMobilityOmega().
std::vector<int> allcells(num_cells); std::vector<int> allcells(num_cells);
for (int cell = 0; cell < num_cells; ++cell) { for (int cell = 0; cell < num_cells; ++cell) {
allcells[cell] = cell; allcells[cell] = cell;
} }
// Warn if any parameters are unused. // Warn if any parameters are unused.
if (param.anyUnused()) { if (param.anyUnused()) {
std::cout << "-------------------- Unused parameters: --------------------\n"; std::cout << "-------------------- Unused parameters: --------------------\n";
param.displayUsage(); param.displayUsage();
std::cout << "----------------------------------------------------------------" << std::endl; std::cout << "----------------------------------------------------------------" << std::endl;
} }
// Write parameters used for later reference. // Write parameters used for later reference.
if (output) { if (output) {
param.writeParam(output_dir + "/spu_2p.param"); param.writeParam(output_dir + "/spu_2p.param");
} }
// Main simulation loop. // Main simulation loop.
@ -450,66 +631,135 @@ main(int argc, char** argv)
Opm::time::StopWatch total_timer; Opm::time::StopWatch total_timer;
total_timer.start(); total_timer.start();
std::cout << "\n\n================ Starting main simulation loop ===============" << std::endl; std::cout << "\n\n================ Starting main simulation loop ===============" << std::endl;
for (int pstep = 0; pstep < num_psteps; ++pstep) { double init_satvol[2] = { 0.0 };
std::cout << "\n\n--------------- Simulation step number " << pstep double init_polymass = 0.0;
<< " ---------------" double satvol[2] = { 0.0 };
<< "\n Current time (days) " << Opm::unit::convert::to(current_time, Opm::unit::day) double polymass = 0.0;
<< "\n Current stepsize (days) " << Opm::unit::convert::to(stepsize, Opm::unit::day) double polymass_adsorbed = 0.0;
<< "\n Total time (days) " << Opm::unit::convert::to(total_time, Opm::unit::day) double injected[2] = { 0.0 };
<< "\n" << std::endl; double produced[2] = { 0.0 };
double polyinj = 0.0;
double polyprod = 0.0;
double tot_injected[2] = { 0.0 };
double tot_produced[2] = { 0.0 };
double tot_polyinj = 0.0;
double tot_polyprod = 0.0;
Opm::computeSaturatedVol(porevol, state.saturation(), init_satvol);
std::cout << "\nInitial saturations are " << init_satvol[0]/tot_porevol
<< " " << init_satvol[1]/tot_porevol << std::endl;
Opm::Watercut watercut;
watercut.push(0.0, 0.0, 0.0);
for (; !simtimer.done(); ++simtimer) {
// Report timestep and (optionally) write state to disk.
simtimer.report(std::cout);
if (output && (simtimer.currentStepNum() % output_interval == 0)) {
outputState(*grid->c_grid(), state, simtimer.currentStepNum(), output_dir);
}
if (output) { // Solve pressure.
outputState(grid->c_grid(), state, pstep, output_dir); if (use_gravity) {
} computeTotalMobilityOmega(*props, polydata, allcells, state.saturation(), state.concentration(),
totmob, omega);
} else {
computeTotalMobility(*props, polydata, allcells, state.saturation(), state.concentration(),
totmob);
}
pressure_timer.start();
psolver.solve(totmob, omega, src, bcs.c_bcs(), state.pressure(), state.faceflux());
pressure_timer.stop();
double pt = pressure_timer.secsSinceStart();
std::cout << "Pressure solver took: " << pt << " seconds." << std::endl;
ptime += pt;
if (use_gravity) { // Process transport sources (to include bdy terms).
computeTotalMobilityOmega(*props, polydata, allcells, state.saturation(), state.concentration(), Opm::computeTransportSource(*grid->c_grid(), src, state.faceflux(), 1.0, reorder_src);
totmob, omega);
} else {
computeTotalMobility(*props, polydata, allcells, state.saturation(), state.concentration(),
totmob);
}
pressure_timer.start();
psolver.solve(totmob, omega, src, bcs.c_bcs(), state.pressure(), state.faceflux());
pressure_timer.stop();
double pt = pressure_timer.secsSinceStart();
std::cout << "Pressure solver took: " << pt << " seconds." << std::endl;
ptime += pt;
const double inflowc0 = poly_inflow(current_time + 1e-5*stepsize); // Find inflow rate.
const double inflowc1 = poly_inflow(current_time + (1.0 - 1e-5)*stepsize); const double current_time = simtimer.currentTime();
if (inflowc0 != inflowc1) { const double stepsize = simtimer.currentStepLength();
std::cout << "**** Warning: polymer inflow rate changes during timestep. Using rate near start of step."; const double inflowc0 = poly_inflow(current_time + 1e-5*stepsize);
} const double inflowc1 = poly_inflow(current_time + (1.0 - 1e-5)*stepsize);
const double inflow_c = inflowc0; if (inflowc0 != inflowc1) {
Opm::toWaterSat(state.saturation(), reorder_sat); std::cout << "**** Warning: polymer inflow rate changes during timestep. Using rate near start of step.";
// We must treat reorder_src here, }
// if we are to handle anything but simple water const double inflow_c = inflowc0;
// 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.
transport_timer.start();
tmodel.solve(&state.faceflux()[0], &reorder_src[0], stepsize, inflow_c,
&reorder_sat[0], &state.concentration()[0], &state.cmax()[0]);
transport_timer.stop();
double tt = transport_timer.secsSinceStart();
std::cout << "Transport solver took: " << tt << " seconds." << std::endl;
ttime += tt;
Opm::toBothSat(reorder_sat, state.saturation());
current_time += stepsize; // Solve transport.
transport_timer.start();
Opm::toWaterSat(state.saturation(), reorder_sat);
tmodel.solve(&state.faceflux()[0], &reorder_src[0], stepsize, inflow_c,
&reorder_sat[0], &state.concentration()[0], &state.cmax()[0]);
Opm::toBothSat(reorder_sat, state.saturation());
transport_timer.stop();
double tt = transport_timer.secsSinceStart();
std::cout << "Transport solver took: " << tt << " seconds." << std::endl;
ttime += tt;
// Report volume balances.
Opm::computeSaturatedVol(porevol, state.saturation(), satvol);
polymass = Opm::computePolymerMass(porevol, state.saturation(), state.concentration(), polydata.deadPoreVol());
polymass_adsorbed = Opm::computePolymerAdsorbed(polydata, porevol, state.cmax());
Opm::computeInjectedProduced(*props, polydata, state.saturation(), state.concentration(),
src, simtimer.currentStepLength(), inflow_c,
injected, produced, polyinj, polyprod);
tot_injected[0] += injected[0];
tot_injected[1] += injected[1];
tot_produced[0] += produced[0];
tot_produced[1] += produced[1];
tot_polyinj += polyinj;
tot_polyprod += polyprod;
std::cout.precision(5);
const int width = 18;
std::cout << "\nVolume and polymer mass balance: "
" water(pv) oil(pv) polymer(kg)\n";
std::cout << " Saturated volumes: "
<< std::setw(width) << satvol[0]/tot_porevol
<< std::setw(width) << satvol[1]/tot_porevol
<< std::setw(width) << polymass << std::endl;
std::cout << " Adsorbed volumes: "
<< std::setw(width) << 0.0
<< std::setw(width) << 0.0
<< std::setw(width) << polymass_adsorbed << std::endl;
std::cout << " Injected volumes: "
<< std::setw(width) << injected[0]/tot_porevol
<< std::setw(width) << injected[1]/tot_porevol
<< std::setw(width) << polyinj << std::endl;
std::cout << " Produced volumes: "
<< std::setw(width) << produced[0]/tot_porevol
<< std::setw(width) << produced[1]/tot_porevol
<< std::setw(width) << polyprod << std::endl;
std::cout << " Total inj volumes: "
<< std::setw(width) << tot_injected[0]/tot_porevol
<< std::setw(width) << tot_injected[1]/tot_porevol
<< std::setw(width) << tot_polyinj << std::endl;
std::cout << " Total prod volumes: "
<< std::setw(width) << tot_produced[0]/tot_porevol
<< std::setw(width) << tot_produced[1]/tot_porevol
<< std::setw(width) << tot_polyprod << 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::setw(width) << (polymass + tot_polyprod - tot_polyinj + polymass_adsorbed) << std::endl;
std::cout << " Init - now - pr + inj: "
<< std::setw(width) << (init_satvol[0] - satvol[0] - tot_produced[0] + tot_injected[0])/tot_porevol
<< std::setw(width) << (init_satvol[1] - satvol[1] - tot_produced[1] + tot_injected[1])/tot_porevol
<< std::setw(width) << (init_polymass - polymass - tot_polyprod + tot_polyinj - polymass_adsorbed)
<< std::endl;
std::cout.precision(8);
watercut.push(simtimer.currentTime() + simtimer.currentStepLength(),
produced[0]/(produced[0] + produced[1]),
tot_produced[0]/tot_porevol);
} }
total_timer.stop(); total_timer.stop();
std::cout << "\n\n================ End of simulation ===============\n" std::cout << "\n\n================ End of simulation ===============\n"
<< "Total time taken: " << total_timer.secsSinceStart() << "Total time taken: " << total_timer.secsSinceStart()
<< "\n Pressure time: " << ptime << "\n Pressure time: " << ptime
<< "\n Transport time: " << ttime << std::endl; << "\n Transport time: " << ttime << std::endl;
if (output) { if (output) {
outputState(grid->c_grid(), state, num_psteps, output_dir); outputState(*grid->c_grid(), state, simtimer.currentStepNum(), output_dir);
outputWaterCut(watercut, output_dir);
} }
} }

110
opm/polymer/polymerUtilities.cpp
View File

@ -98,5 +98,115 @@ namespace Opm
} }
/// @brief Computes injected and produced volumes of all phases,
/// and injeced and produced polymer mass.
/// 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] polyprops polymer properties
/// @param[in] s saturation values (for all P phases)
/// @param[in] c polymer concentration
/// @param[in] src if < 0: total outflow, if > 0: first phase inflow.
/// @param[in] dt timestep used
/// @param[in] inj_c injected concentration
/// @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.
/// @param[out] polyinj injected mass of polymer
/// @param[out] polyprod produced mass of polymer
void computeInjectedProduced(const IncompPropertiesInterface& props,
const Opm::PolymerProperties& polyprops,
const std::vector<double>& s,
const std::vector<double>& c,
const std::vector<double>& src,
const double dt,
const double inj_c,
double* injected,
double* produced,
double& polyinj,
double& polyprod)
{
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);
polyinj = 0.0;
polyprod = 0.0;
std::vector<double> inv_eff_visc(np);
const double* visc = props.viscosity();
std::vector<double> mob(np);
for (int cell = 0; cell < num_cells; ++cell) {
if (src[cell] > 0.0) {
injected[0] += src[cell]*dt;
polyinj += src[cell]*dt*inj_c;
} else if (src[cell] < 0.0) {
const double flux = -src[cell]*dt;
const double* sat = &s[np*cell];
props.relperm(1, sat, &cell, &mob[0], 0);
polyprops.effectiveInvVisc(c[cell], visc, &inv_eff_visc[0]);
double totmob = 0.0;
for (int p = 0; p < np; ++p) {
mob[p] *= inv_eff_visc[p];
totmob += mob[p];
}
for (int p = 0; p < np; ++p) {
produced[p] += (mob[p]/totmob)*flux;
}
polyprod += (mob[0]/totmob)*flux*c[cell]; // TODO check this term.
}
}
}
/// @brief Computes total polymer mass over all grid cells.
/// @param[in] pv the pore volume by cell.
/// @param[in] s saturation values (for all P phases)
/// @param[in] c polymer concentration
/// @param[in] dps dead pore space
/// @return total polymer mass in grid.
double computePolymerMass(const std::vector<double>& pv,
const std::vector<double>& s,
const std::vector<double>& c,
const double dps)
{
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.");
}
double polymass = 0.0;
for (int cell = 0; cell < num_cells; ++cell) {
polymass += c[cell]*pv[cell]*(s[np*cell + 0] - dps);
}
return polymass;
}
/// @brief Computes total absorbed polymer mass over all grid cells.
/// @param[in] polyprops polymer properties
/// @param[in] pv the pore volume by cell.
/// @param[in] cmax max polymer concentration for cell
/// @return total absorbed polymer mass.
double computePolymerAdsorbed(const Opm::PolymerProperties& polyprops,
const std::vector<double>& pv,
const std::vector<double>& cmax)
{
const int num_cells = pv.size();
const double rhor = polyprops.rockDensity();
double abs_mass = 0.0;
for (int cell = 0; cell < num_cells; ++cell) {
abs_mass += polyprops.adsorbtion(cmax[cell])*pv[cell]*rhor;
}
return abs_mass;
}
} // namespace Opm } // namespace Opm

52
opm/polymer/polymerUtilities.hpp
View File

@ -60,6 +60,58 @@ namespace Opm
const std::vector<double>& c, const std::vector<double>& c,
std::vector<double>& totmob, std::vector<double>& totmob,
std::vector<double>& omega); std::vector<double>& omega);
/// @brief Computes injected and produced volumes of all phases,
/// and injeced and produced polymer mass.
/// 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] polyprops polymer properties
/// @param[in] s saturation values (for all P phases)
/// @param[in] c polymer concentration
/// @param[in] src if < 0: total outflow, if > 0: first phase inflow.
/// @param[in] dt timestep used
/// @param[in] inj_c injected concentration
/// @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.
/// @param[out] polyinj injected mass of polymer
/// @param[out] polyprod produced mass of polymer
void computeInjectedProduced(const IncompPropertiesInterface& props,
const Opm::PolymerProperties& polyprops,
const std::vector<double>& s,
const std::vector<double>& c,
const std::vector<double>& src,
const double dt,
const double inj_c,
double* injected,
double* produced,
double& polyinj,
double& polyprod);
/// @brief Computes total (free) polymer mass over all grid cells.
/// @param[in] pv the pore volume by cell.
/// @param[in] s saturation values (for all P phases)
/// @param[in] c polymer concentration
/// @param[in] dps dead pore space
/// @return total polymer mass in grid.
double computePolymerMass(const std::vector<double>& pv,
const std::vector<double>& s,
const std::vector<double>& c,
const double dps);
/// @brief Computes total absorbed polymer mass over all grid cells.
/// @param[in] polyprops polymer properties
/// @param[in] pv the pore volume by cell.
/// @param[in] cmax max polymer concentration for cell
/// @return total absorbed polymer mass.
double computePolymerAdsorbed(const Opm::PolymerProperties& polyprops,
const std::vector<double>& pv,
const std::vector<double>& cmax);
} // namespace Opm } // namespace Opm