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

Make polymer inflow specification by cell.

Browse Source 
Before, it was a single number for the whole domain. It could vary
with time, though.

Simulator behaviour should be identical before and after this commit.

This touches the classes TransportModel*Polymer, Simulator*Polymer,
and the computeInjectedProduced() functions. The PolymerInflow class
and its usage has been replaced with PolymerInflowInterface, and
a subclass PolymerInflowBasic has been created which provides the
old behaviour (using parameters poly_start_days, poly_end_days and
poly_amount).
This commit is contained in:
Atgeirr Flø Rasmussen 2012-10-04 13:38:22 +02:00
parent f5eb116410
commit 6033e5ca69
9 changed files with 91 additions and 71 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

3
examples/test_singlecellsolves.cpp
View File

@ -211,6 +211,7 @@ main(int argc, char** argv)
const double ff = s; // Simplified a lot...
for (int conc = 0; conc < num_concs; ++conc) {
const double c = poly_props.cMax()*double(conc)/double(num_concs - 1);
std::vector<double> polymer_inflow_c(num_cells, c);
// std::cout << "(s, c) = (" << s << ", " << c << ")\n";
transport_src[0] = src[0]*ff;
// Resetting the state for next run.
@ -223,8 +224,8 @@ main(int argc, char** argv)
reorder_model.solve(&state.faceflux()[0],
&porevol[0],
&transport_src[0],
&polymer_inflow_c[0],
dt,
c,
state.saturation(),
state.concentration(),
state.maxconcentration());

26
opm/polymer/SimulatorCompressiblePolymer.cpp
View File

@ -129,7 +129,7 @@ namespace Opm
std::vector< std::vector<int> > columns_;
// Misc. data
std::vector<int> allcells_;
PolymerInflow poly_inflow_;
boost::scoped_ptr<PolymerInflowInterface> poly_inflow_;
};
@ -191,10 +191,7 @@ namespace Opm
tsolver_(grid, props, poly_props, *rock_comp_props,
TransportModelCompressiblePolymer::Bracketing,
param.getDefault("nl_tolerance", 1e-9),
param.getDefault("nl_maxiter", 30)),
poly_inflow_(param.getDefault("poly_start_days", 300.0)*Opm::unit::day,
param.getDefault("poly_end_days", 800.0)*Opm::unit::day,
param.getDefault("poly_amount", poly_props.cMax()))
param.getDefault("nl_maxiter", 30))
{
// For output.
output_ = param.getDefault("output", true);
@ -230,6 +227,11 @@ namespace Opm
extractColumn(grid_, columns_);
}
// Polymer inflow control.
poly_inflow_.reset(new PolymerInflowBasic(param.getDefault("poly_start_days", 300.0)*Opm::unit::day,
param.getDefault("poly_end_days", 800.0)*Opm::unit::day,
param.getDefault("poly_amount", poly_props.cMax())));
// Misc init.
const int num_cells = grid.number_of_cells;
allcells_.resize(num_cells);
@ -244,7 +246,8 @@ namespace Opm
PolymerBlackoilState& state,
WellState& well_state)
{
std::vector<double> transport_src;
std::vector<double> transport_src(grid_.number_of_cells);
std::vector<double> polymer_inflow_c(grid_.number_of_cells);
// Initialisation.
std::vector<double> initial_pressure;
@ -320,12 +323,7 @@ namespace Opm
// Find inflow rate.
const double current_time = timer.currentTime();
double stepsize = timer.currentStepLength();
const double inflowc0 = poly_inflow_(current_time + 1e-5*stepsize);
const double inflowc1 = poly_inflow_(current_time + (1.0 - 1e-5)*stepsize);
if (inflowc0 != inflowc1) {
std::cout << "**** Warning: polymer inflow rate changes during timestep. Using rate near start of step.";
}
const double inflow_c = inflowc0;
poly_inflow_->getInflowValues(current_time, current_time + stepsize, polymer_inflow_c);
// Solve transport.
transport_timer.start();
@ -340,7 +338,7 @@ namespace Opm
for (int tr_substep = 0; tr_substep < num_transport_substeps_; ++tr_substep) {
tsolver_.solve(&state.faceflux()[0], initial_pressure,
state.pressure(), &initial_porevol[0], &porevol[0],
&transport_src[0], stepsize, inflow_c,
&transport_src[0], &polymer_inflow_c[0], stepsize,
state.saturation(), state.surfacevol(),
state.concentration(), state.maxconcentration());
double substep_injected[2] = { 0.0 };
@ -350,7 +348,7 @@ namespace Opm
Opm::computeInjectedProduced(props_, poly_props_,
state.pressure(), state.surfacevol(), state.saturation(),
state.concentration(), state.maxconcentration(),
transport_src, stepsize, inflow_c,
transport_src, polymer_inflow_c, stepsize,
substep_injected, substep_produced,
substep_polyinj, substep_polyprod);
injected[0] += substep_injected[0];

26
opm/polymer/SimulatorPolymer.cpp
View File

@ -126,7 +126,7 @@ namespace Opm
std::vector< std::vector<int> > columns_;
// Misc. data
std::vector<int> allcells_;
PolymerInflow poly_inflow_;
boost::scoped_ptr<PolymerInflowInterface> poly_inflow_;
};
@ -187,10 +187,7 @@ namespace Opm
gravity, wells, src, bcs),
tsolver_(grid, props, poly_props, TransportModelPolymer::Bracketing,
param.getDefault("nl_tolerance", 1e-9),
param.getDefault("nl_maxiter", 30)),
poly_inflow_(param.getDefault("poly_start_days", 300.0)*Opm::unit::day,
param.getDefault("poly_end_days", 800.0)*Opm::unit::day,
param.getDefault("poly_amount", poly_props.cMax()))
param.getDefault("nl_maxiter", 30))
{
// For output.
output_ = param.getDefault("output", true);
@ -226,6 +223,11 @@ namespace Opm
extractColumn(grid_, columns_);
}
// Polymer inflow control.
poly_inflow_.reset(new PolymerInflowBasic(param.getDefault("poly_start_days", 300.0)*Opm::unit::day,
param.getDefault("poly_end_days", 800.0)*Opm::unit::day,
param.getDefault("poly_amount", poly_props.cMax())));
// Misc init.
const int num_cells = grid.number_of_cells;
allcells_.resize(num_cells);
@ -240,7 +242,8 @@ namespace Opm
PolymerState& state,
WellState& well_state)
{
std::vector<double> transport_src;
std::vector<double> transport_src(grid_.number_of_cells);
std::vector<double> polymer_inflow_c(grid_.number_of_cells);
// Initialisation.
std::vector<double> porevol;
@ -316,12 +319,7 @@ namespace Opm
// Find inflow rate.
const double current_time = timer.currentTime();
double stepsize = timer.currentStepLength();
const double inflowc0 = poly_inflow_(current_time + 1e-5*stepsize);
const double inflowc1 = poly_inflow_(current_time + (1.0 - 1e-5)*stepsize);
if (inflowc0 != inflowc1) {
std::cout << "**** Warning: polymer inflow rate changes during timestep. Using rate near start of step.";
}
const double inflow_c = inflowc0;
poly_inflow_->getInflowValues(current_time, current_time + stepsize, polymer_inflow_c);
// Solve transport.
transport_timer.start();
@ -335,11 +333,11 @@ namespace Opm
double substep_polyprod = 0.0;
injected[0] = injected[1] = produced[0] = produced[1] = polyinj = polyprod = 0.0;
for (int tr_substep = 0; tr_substep < num_transport_substeps_; ++tr_substep) {
tsolver_.solve(&state.faceflux()[0], &porevol[0], &transport_src[0], stepsize, inflow_c,
tsolver_.solve(&state.faceflux()[0], &porevol[0], &transport_src[0], &polymer_inflow_c[0], stepsize,
state.saturation(), state.concentration(), state.maxconcentration());
Opm::computeInjectedProduced(props_, poly_props_,
state.saturation(), state.concentration(), state.maxconcentration(),
transport_src, stepsize, inflow_c,
transport_src, polymer_inflow_c, stepsize,
substep_injected, substep_produced, substep_polyinj, substep_polyprod);
injected[0] += substep_injected[0];
injected[1] += substep_injected[1];

10
opm/polymer/TransportModelCompressiblePolymer.cpp
View File

@ -168,8 +168,8 @@ namespace Opm
porevolume0_(0),
porevolume_(0),
source_(0),
polymer_inflow_c_(0),
dt_(0.0),
inflow_c_(0.0),
tol_(tol),
maxit_(maxit),
method_(method),
@ -220,8 +220,8 @@ namespace Opm
const double* porevolume0,
const double* porevolume,
const double* source,
const double* polymer_inflow_c,
const double dt,
const double inflow_c,
std::vector<double>& saturation,
std::vector<double>& surfacevol,
std::vector<double>& concentration,
@ -232,7 +232,7 @@ namespace Opm
porevolume_ = porevolume;
source_ = source;
dt_ = dt;
inflow_c_ = inflow_c;
polymer_inflow_c_ = polymer_inflow_c;
toWaterSat(saturation, saturation_);
concentration_ = &concentration[0];
cmax_ = &cmax[0];
@ -363,7 +363,7 @@ namespace Opm
bool src_is_inflow = src_flux < 0.0;
B_cell0 = 1.0/tm.A0_[np*np*cell + 0];
B_cell = 1.0/tm.A_[np*np*cell + 0];
// Not clear why we multiply by B_cell source terms.
// influx = src_is_inflow ? B_cell*src_flux : 0.0; // Use this after changing transport source.
influx = src_is_inflow ? src_flux : 0.0;
outflux = !src_is_inflow ? src_flux : 0.0;
porevolume0 = tm.porevolume0_[cell];
@ -376,7 +376,7 @@ namespace Opm
rhor = tm.polyprops_.rockDensity();
tm.polyprops_.adsorption(c0, cmax0, ads0);
double mc;
tm.computeMc(tm.inflow_c_, mc);
tm.computeMc(tm.polymer_inflow_c_[cell_index], mc);
influx_polymer = src_is_inflow ? src_flux*mc : 0.0;
for (int i = tm.grid_.cell_facepos[cell]; i < tm.grid_.cell_facepos[cell+1]; ++i) {
int f = tm.grid_.cell_faces[i];

13
opm/polymer/TransportModelCompressiblePolymer.hpp
View File

@ -78,9 +78,14 @@ namespace Opm
/// \param[in] pressure Array with pressure.
/// \param[in] porevolume0 Array with pore volume at start of timestep.
/// \param[in] porevolume Array with pore volume.
/// \param[in] source Transport source term.
/// \param[in] source Transport source term, to be interpreted by sign:
/// (+) Inflow, value is first phase flow (water)
/// per second, in *surface* volumes (unlike the
/// incompressible version).
/// (-) Outflow, value is total flow of all phases
/// per second, in reservoir volumes.
/// \param[in] polymer_inflow_c Array of inflow polymer concentrations per cell.
/// \param[in] dt Time step.
/// \param[in] inflow_c Inflow polymer.
/// \param[in, out] saturation Phase saturations.
/// \param[in, out] surfacevol Surface volumes.
/// \param[in, out] concentration Polymer concentration.
@ -91,8 +96,8 @@ namespace Opm
const double* porevolume0,
const double* porevolume,
const double* source,
const double* polymer_inflow_c,
const double dt,
const double inflow_c,
std::vector<double>& saturation,
std::vector<double>& surfacevol,
std::vector<double>& concentration,
@ -134,8 +139,8 @@ namespace Opm
const double* porevolume0_; // one volume per cell
const double* porevolume_; // one volume per cell
const double* source_; // one source per cell
const double* polymer_inflow_c_;
double dt_;
double inflow_c_;
double tol_;
double maxit_;
SingleCellMethod method_;

8
opm/polymer/TransportModelPolymer.cpp
View File

@ -189,8 +189,8 @@ namespace Opm
maxit_(maxit),
darcyflux_(0),
source_(0),
polymer_inflow_c_(0),
dt_(0.0),
inflow_c_(0.0),
concentration_(0),
cmax_(0),
fractionalflow_(grid.number_of_cells, -1.0),
@ -232,8 +232,8 @@ namespace Opm
void TransportModelPolymer::solve(const double* darcyflux,
const double* porevolume,
const double* source,
const double* polymer_inflow_c,
const double dt,
const double inflow_c,
std::vector<double>& saturation,
std::vector<double>& concentration,
std::vector<double>& cmax)
@ -241,8 +241,8 @@ namespace Opm
darcyflux_ = darcyflux;
porevolume_ = porevolume;
source_ = source;
polymer_inflow_c_ = polymer_inflow_c;
dt_ = dt;
inflow_c_ = inflow_c;
toWaterSat(saturation, saturation_);
concentration_ = &concentration[0];
cmax_ = &cmax[0];
@ -350,7 +350,7 @@ namespace Opm
bool src_is_inflow = dflux < 0.0;
influx = src_is_inflow ? dflux : 0.0;
double mc;
tm.computeMc(tm.inflow_c_, mc);
tm.computeMc(tm.polymer_inflow_c_[cell_index], mc);
influx_polymer = src_is_inflow ? dflux*mc : 0.0;
outflux = !src_is_inflow ? dflux : 0.0;
comp_term = tm.source_[cell]; // Note: this assumes that all source flux is water.

12
opm/polymer/TransportModelPolymer.hpp
View File

@ -67,17 +67,21 @@ namespace Opm
/// Using implicit Euler scheme, reordered.
/// \param[in] darcyflux Array of signed face fluxes.
/// \param[in] porevolume Array of pore volumes.
/// \param[in] source Transport source term.
/// \param[in] source Transport source term, to be interpreted by sign:
/// (+) Inflow, value is first phase flow (water)
/// per second, in reservoir volumes.
/// (-) Outflow, value is total flow of all phases
/// per second, in reservoir volumes.
/// \param[in] polymer_inflow_c Array of inflow polymer concentrations per cell.
/// \param[in] dt Time step.
/// \param[in] inflow_c Time step.
/// \param[in, out] saturation Phase saturations.
/// \param[in, out] concentration Polymer concentration.
/// \param[in, out] cmax Highest concentration that has occured in a given cell.
void solve(const double* darcyflux,
const double* porevolume,
const double* source,
const double* polymer_inflow_c,
const double dt,
const double inflow_c,
std::vector<double>& saturation,
std::vector<double>& concentration,
std::vector<double>& cmax);
@ -149,8 +153,8 @@ namespace Opm
const double* darcyflux_; // one flux per grid face
const double* source_; // one source per cell
const double* polymer_inflow_c_;
double dt_;
double inflow_c_;
std::vector<double> saturation_; // one per cell, only water saturation!
double* concentration_;
double* cmax_;

14
opm/polymer/polymerUtilities.cpp
View File

@ -103,8 +103,8 @@ namespace Opm
/// @param[in] c polymer concentration
/// @param[in] cmax polymer maximum concentration
/// @param[in] src if < 0: total outflow, if > 0: first phase inflow.
/// @param[in] inj_c injected concentration by cell
/// @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.
@ -116,8 +116,8 @@ namespace Opm
const std::vector<double>& c,
const std::vector<double>& cmax,
const std::vector<double>& src,
const std::vector<double>& inj_c,
const double dt,
const double inj_c,
double* injected,
double* produced,
double& polyinj,
@ -139,7 +139,7 @@ namespace Opm
for (int cell = 0; cell < num_cells; ++cell) {
if (src[cell] > 0.0) {
injected[0] += src[cell]*dt;
polyinj += src[cell]*dt*inj_c;
polyinj += src[cell]*dt*inj_c[cell];
} else if (src[cell] < 0.0) {
const double flux = -src[cell]*dt;
const double* sat = &s[np*cell];
@ -170,15 +170,13 @@ namespace Opm
/// @param[in] c polymer concentration
/// @param[in] cmax polymer maximum concentration
/// @param[in] src if < 0: total outflow, if > 0: first phase inflow.
/// @param[in] inj_c injected concentration by cell
/// @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 BlackoilPropertiesInterface& props,
const Opm::PolymerProperties& polyprops,
const std::vector<double>& press,
@ -187,8 +185,8 @@ namespace Opm
const std::vector<double>& c,
const std::vector<double>& cmax,
const std::vector<double>& src,
const std::vector<double>& inj_c,
const double dt,
const double inj_c,
double* injected,
double* produced,
double& polyinj,
@ -210,7 +208,7 @@ namespace Opm
for (int cell = 0; cell < num_cells; ++cell) {
if (src[cell] > 0.0) {
injected[0] += src[cell]*dt;
polyinj += src[cell]*dt*inj_c;
polyinj += src[cell]*dt*inj_c[cell];
} else if (src[cell] < 0.0) {
const double flux = -src[cell]*dt;
const double* sat = &s[np*cell];

50
opm/polymer/polymerUtilities.hpp
View File

@ -77,8 +77,8 @@ namespace Opm
/// @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] inj_c injected concentration by cell
/// @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.
@ -90,8 +90,8 @@ namespace Opm
const std::vector<double>& c,
const std::vector<double>& cmax,
const std::vector<double>& src,
const std::vector<double>& inj_c,
const double dt,
const double inj_c,
double* injected,
double* produced,
double& polyinj,
@ -111,8 +111,8 @@ namespace Opm
/// @param[in] c polymer concentration
/// @param[in] cmax polymer maximum concentration
/// @param[in] src if < 0: total outflow, if > 0: first phase inflow.
/// @param[in] inj_c injected concentration by cell
/// @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().
@ -127,9 +127,9 @@ namespace Opm
const std::vector<double>& s,
const std::vector<double>& c,
const std::vector<double>& cmax,
const std::vector<double>& src,
const double dt,
const double inj_c,
const std::vector<double>& src,
const std::vector<double>& inj_c,
const double dt,
double* injected,
double* produced,
double& polyinj,
@ -172,28 +172,44 @@ namespace Opm
const RockCompressibility* rock_comp);
/// @brief Functor giving the injected amount of polymer as a function of time.
class PolymerInflow
class PolymerInflowInterface
{
public:
virtual ~PolymerInflowInterface() {}
virtual void getInflowValues(const double step_start,
const double step_end,
std::vector<double>& poly_inflow_c) = 0;
};
/// @brief Functor giving the injected amount of polymer per cell as a function of time.
class PolymerInflowBasic : public PolymerInflowInterface
{
public:
/// Constructor.
/// @param[in] starttime Start time of injection in seconds.
/// @param[in] endtime End time of injection in seconds.
/// @param[in] amount Amount to be injected per second.
PolymerInflow(const double starttime,
const double endtime,
const double amount)
PolymerInflowBasic(const double starttime,
const double endtime,
const double amount)
: stime_(starttime), etime_(endtime), amount_(amount)
{
}
/// Get the current injection rate.
/// @param[in] time Current time in seconds.
double operator()(double time)
virtual void getInflowValues(const double step_start,
const double step_end,
std::vector<double>& poly_inflow_c)
{
if (time >= stime_ && time < etime_) {
return amount_;
const double eps = 1e-5*(step_end - step_start);
if (step_start + eps >= stime_ && step_end - eps <= etime_) {
std::fill(poly_inflow_c.begin(), poly_inflow_c.end(), amount_);
} else if (step_start + eps <= etime_ && step_end - eps >= stime_) {
MESSAGE("Warning: polymer injection set to change inside timestep. Using value at start of step.");
std::fill(poly_inflow_c.begin(), poly_inflow_c.end(), amount_);
} else {
return 0.0;
std::fill(poly_inflow_c.begin(), poly_inflow_c.end(), 0.0);
}
}
private: