mirror of
https://github.com/OPM/opm-simulators.git
synced 2025-02-25 18:55:30 -06:00
Merge remote-tracking branch 'upstream/master'
This commit is contained in:
@@ -28,8 +28,8 @@ namespace Opm
|
|||||||
{
|
{
|
||||||
rock_.init(deck, grid);
|
rock_.init(deck, grid);
|
||||||
pvt_.init(deck, 200);
|
pvt_.init(deck, 200);
|
||||||
SaturationPropsFromDeck<SatFuncStone2Uniform>* ptr
|
SaturationPropsFromDeck<SatFuncSimpleUniform>* ptr
|
||||||
= new SaturationPropsFromDeck<SatFuncStone2Uniform>();
|
= new SaturationPropsFromDeck<SatFuncSimpleUniform>();
|
||||||
satprops_.reset(ptr);
|
satprops_.reset(ptr);
|
||||||
ptr->init(deck, grid, 200);
|
ptr->init(deck, grid, 200);
|
||||||
|
|
||||||
@@ -50,30 +50,43 @@ namespace Opm
|
|||||||
// Unfortunate lack of pointer smartness here...
|
// Unfortunate lack of pointer smartness here...
|
||||||
const int sat_samples = param.getDefault("sat_tab_size", 200);
|
const int sat_samples = param.getDefault("sat_tab_size", 200);
|
||||||
std::string threephase_model = param.getDefault<std::string>("threephase_model", "simple");
|
std::string threephase_model = param.getDefault<std::string>("threephase_model", "simple");
|
||||||
bool use_stone2 = (threephase_model == "stone2");
|
|
||||||
if (sat_samples > 1) {
|
if (sat_samples > 1) {
|
||||||
if (use_stone2) {
|
if (threephase_model == "stone2") {
|
||||||
SaturationPropsFromDeck<SatFuncStone2Uniform>* ptr
|
SaturationPropsFromDeck<SatFuncStone2Uniform>* ptr
|
||||||
= new SaturationPropsFromDeck<SatFuncStone2Uniform>();
|
= new SaturationPropsFromDeck<SatFuncStone2Uniform>();
|
||||||
satprops_.reset(ptr);
|
satprops_.reset(ptr);
|
||||||
ptr->init(deck, grid, sat_samples);
|
ptr->init(deck, grid, sat_samples);
|
||||||
} else {
|
} else if (threephase_model == "simple") {
|
||||||
SaturationPropsFromDeck<SatFuncSimpleUniform>* ptr
|
SaturationPropsFromDeck<SatFuncSimpleUniform>* ptr
|
||||||
= new SaturationPropsFromDeck<SatFuncSimpleUniform>();
|
= new SaturationPropsFromDeck<SatFuncSimpleUniform>();
|
||||||
satprops_.reset(ptr);
|
satprops_.reset(ptr);
|
||||||
ptr->init(deck, grid, sat_samples);
|
ptr->init(deck, grid, sat_samples);
|
||||||
|
} else if (threephase_model == "gwseg") {
|
||||||
|
SaturationPropsFromDeck<SatFuncGwsegUniform>* ptr
|
||||||
|
= new SaturationPropsFromDeck<SatFuncGwsegUniform>();
|
||||||
|
satprops_.reset(ptr);
|
||||||
|
ptr->init(deck, grid, sat_samples);
|
||||||
|
} else {
|
||||||
|
THROW("Unknown threephase_model: " << threephase_model);
|
||||||
}
|
}
|
||||||
} else {
|
} else {
|
||||||
if (use_stone2) {
|
if (threephase_model == "stone2") {
|
||||||
SaturationPropsFromDeck<SatFuncStone2Nonuniform>* ptr
|
SaturationPropsFromDeck<SatFuncStone2Nonuniform>* ptr
|
||||||
= new SaturationPropsFromDeck<SatFuncStone2Nonuniform>();
|
= new SaturationPropsFromDeck<SatFuncStone2Nonuniform>();
|
||||||
satprops_.reset(ptr);
|
satprops_.reset(ptr);
|
||||||
ptr->init(deck, grid, sat_samples);
|
ptr->init(deck, grid, sat_samples);
|
||||||
} else {
|
} else if (threephase_model == "simple") {
|
||||||
SaturationPropsFromDeck<SatFuncSimpleNonuniform>* ptr
|
SaturationPropsFromDeck<SatFuncSimpleNonuniform>* ptr
|
||||||
= new SaturationPropsFromDeck<SatFuncSimpleNonuniform>();
|
= new SaturationPropsFromDeck<SatFuncSimpleNonuniform>();
|
||||||
satprops_.reset(ptr);
|
satprops_.reset(ptr);
|
||||||
ptr->init(deck, grid, sat_samples);
|
ptr->init(deck, grid, sat_samples);
|
||||||
|
} else if (threephase_model == "gwseg") {
|
||||||
|
SaturationPropsFromDeck<SatFuncGwsegNonuniform>* ptr
|
||||||
|
= new SaturationPropsFromDeck<SatFuncGwsegNonuniform>();
|
||||||
|
satprops_.reset(ptr);
|
||||||
|
ptr->init(deck, grid, sat_samples);
|
||||||
|
} else {
|
||||||
|
THROW("Unknown threephase_model: " << threephase_model);
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
|
|
||||||
|
|||||||
@@ -26,6 +26,7 @@
|
|||||||
#include <opm/core/fluid/blackoil/BlackoilPhases.hpp>
|
#include <opm/core/fluid/blackoil/BlackoilPhases.hpp>
|
||||||
#include <opm/core/fluid/SatFuncStone2.hpp>
|
#include <opm/core/fluid/SatFuncStone2.hpp>
|
||||||
#include <opm/core/fluid/SatFuncSimple.hpp>
|
#include <opm/core/fluid/SatFuncSimple.hpp>
|
||||||
|
#include <opm/core/fluid/SatFuncGwseg.hpp>
|
||||||
#include <vector>
|
#include <vector>
|
||||||
|
|
||||||
struct UnstructuredGrid;
|
struct UnstructuredGrid;
|
||||||
|
|||||||
@@ -26,6 +26,10 @@
|
|||||||
|
|
||||||
#include <opm/core/utility/have_boost_redef.hpp>
|
#include <opm/core/utility/have_boost_redef.hpp>
|
||||||
|
|
||||||
|
// Silence compatibility warning from DUNE headers since we don't use
|
||||||
|
// the deprecated member anyway (in this compilation unit)
|
||||||
|
#define DUNE_COMMON_FIELDVECTOR_SIZE_IS_METHOD 1
|
||||||
|
|
||||||
// TODO: clean up includes.
|
// TODO: clean up includes.
|
||||||
#include <dune/common/deprecated.hh>
|
#include <dune/common/deprecated.hh>
|
||||||
#include <dune/istl/bvector.hh>
|
#include <dune/istl/bvector.hh>
|
||||||
|
|||||||
@@ -36,6 +36,7 @@
|
|||||||
#include <cmath>
|
#include <cmath>
|
||||||
#include <iostream>
|
#include <iostream>
|
||||||
#include <iomanip>
|
#include <iomanip>
|
||||||
|
#include <numeric>
|
||||||
|
|
||||||
namespace Opm
|
namespace Opm
|
||||||
{
|
{
|
||||||
@@ -122,7 +123,7 @@ namespace Opm
|
|||||||
WellState& well_state)
|
WellState& well_state)
|
||||||
{
|
{
|
||||||
const int nc = grid_.number_of_cells;
|
const int nc = grid_.number_of_cells;
|
||||||
const int nw = wells_->number_of_wells;
|
const int nw = (wells_ != 0) ? wells_->number_of_wells : 0;
|
||||||
|
|
||||||
// Set up dynamic data.
|
// Set up dynamic data.
|
||||||
computePerSolveDynamicData(dt, state, well_state);
|
computePerSolveDynamicData(dt, state, well_state);
|
||||||
@@ -446,30 +447,54 @@ namespace Opm
|
|||||||
/// Compute per-iteration dynamic properties for wells.
|
/// Compute per-iteration dynamic properties for wells.
|
||||||
void CompressibleTpfa::computeWellDynamicData(const double /*dt*/,
|
void CompressibleTpfa::computeWellDynamicData(const double /*dt*/,
|
||||||
const BlackoilState& /*state*/,
|
const BlackoilState& /*state*/,
|
||||||
const WellState& /*well_state*/)
|
const WellState& well_state)
|
||||||
{
|
{
|
||||||
// These are the variables that get computed by this function:
|
// These are the variables that get computed by this function:
|
||||||
//
|
//
|
||||||
// std::vector<double> wellperf_A_;
|
// std::vector<double> wellperf_A_;
|
||||||
// std::vector<double> wellperf_phasemob_;
|
// std::vector<double> wellperf_phasemob_;
|
||||||
const int np = props_.numPhases();
|
const int np = props_.numPhases();
|
||||||
const int nw = wells_->number_of_wells;
|
const int nw = (wells_ != 0) ? wells_->number_of_wells : 0;
|
||||||
const int nperf = wells_->well_connpos[nw];
|
const int nperf = (wells_ != 0) ? wells_->well_connpos[nw] : 0;
|
||||||
wellperf_A_.resize(nperf*np*np);
|
wellperf_A_.resize(nperf*np*np);
|
||||||
wellperf_phasemob_.resize(nperf*np);
|
wellperf_phasemob_.resize(nperf*np);
|
||||||
// The A matrix is set equal to the perforation grid cells'
|
// The A matrix is set equal to the perforation grid cells'
|
||||||
// matrix, for both injectors and producers.
|
// matrix for producers, computed from bhp and injection
|
||||||
|
// component fractions from
|
||||||
// The mobilities are set equal to the perforation grid cells'
|
// The mobilities are set equal to the perforation grid cells'
|
||||||
// mobilities, for both injectors and producers.
|
// mobilities for producers.
|
||||||
|
std::vector<double> mu(np);
|
||||||
for (int w = 0; w < nw; ++w) {
|
for (int w = 0; w < nw; ++w) {
|
||||||
|
bool producer = (wells_->type[w] == PRODUCER);
|
||||||
|
const double* comp_frac = &wells_->comp_frac[np*w];
|
||||||
for (int j = wells_->well_connpos[w]; j < wells_->well_connpos[w+1]; ++j) {
|
for (int j = wells_->well_connpos[w]; j < wells_->well_connpos[w+1]; ++j) {
|
||||||
const int c = wells_->well_cells[j];
|
const int c = wells_->well_cells[j];
|
||||||
const double* cA = &cell_A_[np*np*c];
|
|
||||||
double* wpA = &wellperf_A_[np*np*j];
|
double* wpA = &wellperf_A_[np*np*j];
|
||||||
std::copy(cA, cA + np*np, wpA);
|
|
||||||
const double* cM = &cell_phasemob_[np*c];
|
|
||||||
double* wpM = &wellperf_phasemob_[np*j];
|
double* wpM = &wellperf_phasemob_[np*j];
|
||||||
std::copy(cM, cM + np, wpM);
|
if (producer) {
|
||||||
|
const double* cA = &cell_A_[np*np*c];
|
||||||
|
std::copy(cA, cA + np*np, wpA);
|
||||||
|
const double* cM = &cell_phasemob_[np*c];
|
||||||
|
std::copy(cM, cM + np, wpM);
|
||||||
|
} else {
|
||||||
|
const double bhp = well_state.bhp()[w];
|
||||||
|
double perf_p = bhp;
|
||||||
|
for (int phase = 0; phase < np; ++phase) {
|
||||||
|
perf_p += wellperf_gpot_[np*j + phase]*comp_frac[phase];
|
||||||
|
}
|
||||||
|
// Hack warning: comp_frac is used as a component
|
||||||
|
// surface-volume variable in calls to matrix() and
|
||||||
|
// viscosity(), but as a saturation in the call to
|
||||||
|
// relperm(). This is probably ok as long as injectors
|
||||||
|
// only inject pure fluids.
|
||||||
|
props_.matrix(1, &perf_p, comp_frac, &c, wpA, NULL);
|
||||||
|
props_.viscosity(1, &perf_p, comp_frac, &c, &mu[0], NULL);
|
||||||
|
ASSERT(std::fabs(std::accumulate(comp_frac, comp_frac + np, 0.0) - 1.0) < 1e-6);
|
||||||
|
props_.relperm (1, comp_frac, &c, wpM , NULL);
|
||||||
|
for (int phase = 0; phase < np; ++phase) {
|
||||||
|
wpM[phase] /= mu[phase];
|
||||||
|
}
|
||||||
|
}
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
@@ -487,9 +512,9 @@ namespace Opm
|
|||||||
const double* z = &state.surfacevol()[0];
|
const double* z = &state.surfacevol()[0];
|
||||||
UnstructuredGrid* gg = const_cast<UnstructuredGrid*>(&grid_);
|
UnstructuredGrid* gg = const_cast<UnstructuredGrid*>(&grid_);
|
||||||
CompletionData completion_data;
|
CompletionData completion_data;
|
||||||
completion_data.gpot = &wellperf_gpot_[0];
|
completion_data.gpot = ! wellperf_gpot_.empty() ? &wellperf_gpot_[0] : 0;
|
||||||
completion_data.A = &wellperf_A_[0];
|
completion_data.A = ! wellperf_A_.empty() ? &wellperf_A_[0] : 0;
|
||||||
completion_data.phasemob = &wellperf_phasemob_[0];
|
completion_data.phasemob = ! wellperf_phasemob_.empty() ? &wellperf_phasemob_[0] : 0;
|
||||||
cfs_tpfa_res_wells wells_tmp;
|
cfs_tpfa_res_wells wells_tmp;
|
||||||
wells_tmp.W = const_cast<Wells*>(wells_);
|
wells_tmp.W = const_cast<Wells*>(wells_);
|
||||||
wells_tmp.data = &completion_data;
|
wells_tmp.data = &completion_data;
|
||||||
@@ -574,9 +599,9 @@ namespace Opm
|
|||||||
{
|
{
|
||||||
UnstructuredGrid* gg = const_cast<UnstructuredGrid*>(&grid_);
|
UnstructuredGrid* gg = const_cast<UnstructuredGrid*>(&grid_);
|
||||||
CompletionData completion_data;
|
CompletionData completion_data;
|
||||||
completion_data.gpot = const_cast<double*>(&wellperf_gpot_[0]);
|
completion_data.gpot = ! wellperf_gpot_.empty() ? const_cast<double*>(&wellperf_gpot_[0]) : 0;
|
||||||
completion_data.A = const_cast<double*>(&wellperf_A_[0]);
|
completion_data.A = ! wellperf_A_.empty() ? const_cast<double*>(&wellperf_A_[0]) : 0;
|
||||||
completion_data.phasemob = const_cast<double*>(&wellperf_phasemob_[0]);
|
completion_data.phasemob = ! wellperf_phasemob_.empty() ? const_cast<double*>(&wellperf_phasemob_[0]) : 0;
|
||||||
cfs_tpfa_res_wells wells_tmp;
|
cfs_tpfa_res_wells wells_tmp;
|
||||||
wells_tmp.W = const_cast<Wells*>(wells_);
|
wells_tmp.W = const_cast<Wells*>(wells_);
|
||||||
wells_tmp.data = &completion_data;
|
wells_tmp.data = &completion_data;
|
||||||
@@ -584,6 +609,9 @@ namespace Opm
|
|||||||
forces.wells = &wells_tmp;
|
forces.wells = &wells_tmp;
|
||||||
forces.src = NULL;
|
forces.src = NULL;
|
||||||
|
|
||||||
|
double* wpress = ! well_state.bhp ().empty() ? & well_state.bhp ()[0] : 0;
|
||||||
|
double* wflux = ! well_state.perfRates().empty() ? & well_state.perfRates()[0] : 0;
|
||||||
|
|
||||||
cfs_tpfa_res_flux(gg,
|
cfs_tpfa_res_flux(gg,
|
||||||
&forces,
|
&forces,
|
||||||
props_.numPhases(),
|
props_.numPhases(),
|
||||||
@@ -592,9 +620,9 @@ namespace Opm
|
|||||||
&face_phasemob_[0],
|
&face_phasemob_[0],
|
||||||
&face_gravcap_[0],
|
&face_gravcap_[0],
|
||||||
&state.pressure()[0],
|
&state.pressure()[0],
|
||||||
&well_state.bhp()[0],
|
wpress,
|
||||||
&state.faceflux()[0],
|
&state.faceflux()[0],
|
||||||
&well_state.perfRates()[0]);
|
wflux);
|
||||||
cfs_tpfa_res_fpress(gg,
|
cfs_tpfa_res_fpress(gg,
|
||||||
props_.numPhases(),
|
props_.numPhases(),
|
||||||
&htrans_[0],
|
&htrans_[0],
|
||||||
@@ -604,6 +632,23 @@ namespace Opm
|
|||||||
&state.pressure()[0],
|
&state.pressure()[0],
|
||||||
&state.faceflux()[0],
|
&state.faceflux()[0],
|
||||||
&state.facepressure()[0]);
|
&state.facepressure()[0]);
|
||||||
|
|
||||||
|
// Compute well perforation pressures (not done by the C code).
|
||||||
|
if (wells_ != 0) {
|
||||||
|
const int nw = wells_->number_of_wells;
|
||||||
|
const int np = props_.numPhases();
|
||||||
|
for (int w = 0; w < nw; ++w) {
|
||||||
|
const double* comp_frac = &wells_->comp_frac[np*w];
|
||||||
|
for (int j = wells_->well_connpos[w]; j < wells_->well_connpos[w+1]; ++j) {
|
||||||
|
const double bhp = well_state.bhp()[w];
|
||||||
|
double perf_p = bhp;
|
||||||
|
for (int phase = 0; phase < np; ++phase) {
|
||||||
|
perf_p += wellperf_gpot_[np*j + phase]*comp_frac[phase];
|
||||||
|
}
|
||||||
|
well_state.perfPress()[j] = perf_p;
|
||||||
|
}
|
||||||
|
}
|
||||||
|
}
|
||||||
}
|
}
|
||||||
|
|
||||||
} // namespace Opm
|
} // namespace Opm
|
||||||
|
|||||||
@@ -142,7 +142,7 @@ impl_allocate(struct UnstructuredGrid *G ,
|
|||||||
nnu = G->number_of_cells;
|
nnu = G->number_of_cells;
|
||||||
nwperf = 0;
|
nwperf = 0;
|
||||||
|
|
||||||
if (wells != NULL) { assert (wells->W != NULL);
|
if ((wells != NULL) && (wells->W != NULL)) {
|
||||||
nnu += wells->W->number_of_wells;
|
nnu += wells->W->number_of_wells;
|
||||||
nwperf = wells->W->well_connpos[ wells->W->number_of_wells ];
|
nwperf = wells->W->well_connpos[ wells->W->number_of_wells ];
|
||||||
}
|
}
|
||||||
@@ -185,13 +185,15 @@ construct_matrix(struct UnstructuredGrid *G ,
|
|||||||
/* ---------------------------------------------------------------------- */
|
/* ---------------------------------------------------------------------- */
|
||||||
{
|
{
|
||||||
int f, c1, c2, w, i, nc, nnu;
|
int f, c1, c2, w, i, nc, nnu;
|
||||||
|
int wells_present;
|
||||||
size_t nnz;
|
size_t nnz;
|
||||||
|
|
||||||
struct CSRMatrix *A;
|
struct CSRMatrix *A;
|
||||||
|
|
||||||
nc = nnu = G->number_of_cells;
|
nc = nnu = G->number_of_cells;
|
||||||
if (wells != NULL) {
|
|
||||||
assert (wells->W != NULL);
|
wells_present = (wells != NULL) && (wells->W != NULL);
|
||||||
|
if (wells_present) {
|
||||||
nnu += wells->W->number_of_wells;
|
nnu += wells->W->number_of_wells;
|
||||||
}
|
}
|
||||||
|
|
||||||
@@ -214,7 +216,7 @@ construct_matrix(struct UnstructuredGrid *G ,
|
|||||||
}
|
}
|
||||||
}
|
}
|
||||||
|
|
||||||
if (wells != NULL) {
|
if (wells_present) {
|
||||||
/* Well <-> cell connections */
|
/* Well <-> cell connections */
|
||||||
struct Wells *W = wells->W;
|
struct Wells *W = wells->W;
|
||||||
|
|
||||||
@@ -252,7 +254,7 @@ construct_matrix(struct UnstructuredGrid *G ,
|
|||||||
}
|
}
|
||||||
}
|
}
|
||||||
|
|
||||||
if (wells != NULL) {
|
if (wells_present) {
|
||||||
/* Fill well <-> cell connections */
|
/* Fill well <-> cell connections */
|
||||||
struct Wells *W = wells->W;
|
struct Wells *W = wells->W;
|
||||||
|
|
||||||
@@ -741,6 +743,29 @@ assemble_completion_to_cell(int c, int wdof, int np, double dt,
|
|||||||
}
|
}
|
||||||
|
|
||||||
|
|
||||||
|
/* ---------------------------------------------------------------------- */
|
||||||
|
static void
|
||||||
|
welleq_coeff_shut(int np, struct cfs_tpfa_res_data *h,
|
||||||
|
double *res, double *w2c, double *w2w)
|
||||||
|
/* ---------------------------------------------------------------------- */
|
||||||
|
{
|
||||||
|
int p;
|
||||||
|
double fwi;
|
||||||
|
const double *dpflux_w;
|
||||||
|
|
||||||
|
/* Sum reservoir phase flux derivatives set by
|
||||||
|
* compute_darcyflux_and_deriv(). */
|
||||||
|
dpflux_w = h->pimpl->flux_work + (1 * np);
|
||||||
|
for (p = 0, fwi = 0.0; p < np; p++) {
|
||||||
|
fwi += dpflux_w[ p ];
|
||||||
|
}
|
||||||
|
|
||||||
|
*res = 0.0;
|
||||||
|
*w2c = 0.0;
|
||||||
|
*w2w = fwi;
|
||||||
|
}
|
||||||
|
|
||||||
|
|
||||||
/* ---------------------------------------------------------------------- */
|
/* ---------------------------------------------------------------------- */
|
||||||
static void
|
static void
|
||||||
welleq_coeff_bhp(int np, double dp, struct cfs_tpfa_res_data *h,
|
welleq_coeff_bhp(int np, double dp, struct cfs_tpfa_res_data *h,
|
||||||
@@ -795,7 +820,34 @@ welleq_coeff_resv(int np, struct cfs_tpfa_res_data *h,
|
|||||||
|
|
||||||
/* ---------------------------------------------------------------------- */
|
/* ---------------------------------------------------------------------- */
|
||||||
static void
|
static void
|
||||||
assemble_completion_to_well(int w, int c, int nc, int np,
|
welleq_coeff_surfrate(int i, int np, struct cfs_tpfa_res_data *h,
|
||||||
|
struct WellControls *ctrl,
|
||||||
|
double *res, double *w2c, double *w2w)
|
||||||
|
/* ---------------------------------------------------------------------- */
|
||||||
|
{
|
||||||
|
int p;
|
||||||
|
double f;
|
||||||
|
const double *pflux, *dpflux_w, *dpflux_c, *distr;
|
||||||
|
|
||||||
|
pflux = h->pimpl->compflux_p + (i * (1 * np));
|
||||||
|
dpflux_w = h->pimpl->compflux_deriv_p + (i * (2 * np));
|
||||||
|
dpflux_c = dpflux_w + (1 * (1 * np));
|
||||||
|
distr = ctrl->distr + (ctrl->current * (1 * np));
|
||||||
|
|
||||||
|
*res = *w2c = *w2w = 0.0;
|
||||||
|
for (p = 0; p < np; p++) {
|
||||||
|
f = distr[ p ];
|
||||||
|
|
||||||
|
*res += f * pflux [ p ];
|
||||||
|
*w2w += f * dpflux_w[ p ];
|
||||||
|
*w2c += f * dpflux_c[ p ];
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
|
||||||
|
/* ---------------------------------------------------------------------- */
|
||||||
|
static void
|
||||||
|
assemble_completion_to_well(int i, int w, int c, int nc, int np,
|
||||||
double pw, double dt,
|
double pw, double dt,
|
||||||
struct cfs_tpfa_res_wells *wells,
|
struct cfs_tpfa_res_wells *wells,
|
||||||
struct cfs_tpfa_res_data *h )
|
struct cfs_tpfa_res_data *h )
|
||||||
@@ -815,23 +867,25 @@ assemble_completion_to_well(int w, int c, int nc, int np,
|
|||||||
ctrl = W->ctrls[ w ];
|
ctrl = W->ctrls[ w ];
|
||||||
|
|
||||||
if (ctrl->current < 0) {
|
if (ctrl->current < 0) {
|
||||||
assert (0); /* Shut wells currently not supported */
|
/* Interpreting a negative current control index to mean a shut well */
|
||||||
|
welleq_coeff_shut(np, h, &res, &w2c, &w2w);
|
||||||
}
|
}
|
||||||
|
else {
|
||||||
|
switch (ctrl->type[ ctrl->current ]) {
|
||||||
|
case BHP :
|
||||||
|
welleq_coeff_bhp(np, pw - ctrl->target[ ctrl->current ],
|
||||||
|
h, &res, &w2c, &w2w);
|
||||||
|
break;
|
||||||
|
|
||||||
switch (ctrl->type[ ctrl->current ]) {
|
case RESERVOIR_RATE:
|
||||||
case BHP :
|
assert (W->number_of_phases == np);
|
||||||
welleq_coeff_bhp(np, pw - ctrl->target[ ctrl->current ],
|
welleq_coeff_resv(np, h, ctrl, &res, &w2c, &w2w);
|
||||||
h, &res, &w2c, &w2w);
|
break;
|
||||||
break;
|
|
||||||
|
|
||||||
case RESERVOIR_RATE:
|
case SURFACE_RATE:
|
||||||
assert (W->number_of_phases == np);
|
welleq_coeff_surfrate(i, np, h, ctrl, &res, &w2c, &w2w);
|
||||||
welleq_coeff_resv(np, h, ctrl, &res, &w2c, &w2w);
|
break;
|
||||||
break;
|
}
|
||||||
|
|
||||||
case SURFACE_RATE:
|
|
||||||
assert (0); /* Surface rate currently not supported */
|
|
||||||
break;
|
|
||||||
}
|
}
|
||||||
|
|
||||||
/* Assemble completion contributions */
|
/* Assemble completion contributions */
|
||||||
@@ -854,7 +908,7 @@ assemble_well_contrib(struct cfs_tpfa_res_wells *wells ,
|
|||||||
struct cfs_tpfa_res_data *h )
|
struct cfs_tpfa_res_data *h )
|
||||||
{
|
{
|
||||||
int w, i, c, np, np2, nc;
|
int w, i, c, np, np2, nc;
|
||||||
int is_neumann;
|
int is_neumann, is_open;
|
||||||
double pw, dp;
|
double pw, dp;
|
||||||
double *WI, *gpot, *pmobp;
|
double *WI, *gpot, *pmobp;
|
||||||
const double *Ac, *dAc;
|
const double *Ac, *dAc;
|
||||||
@@ -876,6 +930,7 @@ assemble_well_contrib(struct cfs_tpfa_res_wells *wells ,
|
|||||||
|
|
||||||
for (w = i = 0; w < W->number_of_wells; w++) {
|
for (w = i = 0; w < W->number_of_wells; w++) {
|
||||||
pw = wpress[ w ];
|
pw = wpress[ w ];
|
||||||
|
is_open = W->ctrls[w]->current >= 0;
|
||||||
|
|
||||||
for (; i < W->well_connpos[w + 1];
|
for (; i < W->well_connpos[w + 1];
|
||||||
i++, gpot += np, pmobp += np) {
|
i++, gpot += np, pmobp += np) {
|
||||||
@@ -888,14 +943,16 @@ assemble_well_contrib(struct cfs_tpfa_res_wells *wells ,
|
|||||||
|
|
||||||
init_completion_contrib(i, np, Ac, dAc, h->pimpl);
|
init_completion_contrib(i, np, Ac, dAc, h->pimpl);
|
||||||
|
|
||||||
assemble_completion_to_cell(c, nc + w, np, dt, h);
|
if (is_open) {
|
||||||
|
assemble_completion_to_cell(c, nc + w, np, dt, h);
|
||||||
|
}
|
||||||
|
|
||||||
/* Prepare for RESV controls */
|
/* Prepare for RESV controls */
|
||||||
compute_darcyflux_and_deriv(np, WI[i], dp, pmobp, gpot,
|
compute_darcyflux_and_deriv(np, WI[i], dp, pmobp, gpot,
|
||||||
h->pimpl->flux_work,
|
h->pimpl->flux_work,
|
||||||
h->pimpl->flux_work + np);
|
h->pimpl->flux_work + np);
|
||||||
|
|
||||||
assemble_completion_to_well(w, c, nc, np, pw, dt, wells, h);
|
assemble_completion_to_well(i, w, c, nc, np, pw, dt, wells, h);
|
||||||
}
|
}
|
||||||
|
|
||||||
ctrl = W->ctrls[ w ];
|
ctrl = W->ctrls[ w ];
|
||||||
@@ -1127,8 +1184,7 @@ cfs_tpfa_res_construct(struct UnstructuredGrid *G ,
|
|||||||
nf = G->number_of_faces;
|
nf = G->number_of_faces;
|
||||||
nwperf = 0;
|
nwperf = 0;
|
||||||
|
|
||||||
if (wells != NULL) {
|
if ((wells != NULL) && (wells->W != NULL)) {
|
||||||
assert (wells->W != NULL);
|
|
||||||
nwperf = wells->W->well_connpos[ wells->W->number_of_wells ];
|
nwperf = wells->W->well_connpos[ wells->W->number_of_wells ];
|
||||||
}
|
}
|
||||||
|
|
||||||
@@ -1194,7 +1250,9 @@ cfs_tpfa_res_assemble(struct UnstructuredGrid *G ,
|
|||||||
assemble_cell_contrib(G, c, h);
|
assemble_cell_contrib(G, c, h);
|
||||||
}
|
}
|
||||||
|
|
||||||
if ((forces != NULL) && (forces->wells != NULL)) {
|
if ((forces != NULL) &&
|
||||||
|
(forces->wells != NULL) &&
|
||||||
|
(forces->wells->W != NULL)) {
|
||||||
compute_well_compflux_and_deriv(forces->wells, cq->nphases,
|
compute_well_compflux_and_deriv(forces->wells, cq->nphases,
|
||||||
cpress, wpress, h->pimpl);
|
cpress, wpress, h->pimpl);
|
||||||
|
|
||||||
@@ -1297,8 +1355,9 @@ cfs_tpfa_res_flux(struct UnstructuredGrid *G ,
|
|||||||
{
|
{
|
||||||
compute_flux(G, np, trans, pmobf, gravcap_f, cpress, fflux);
|
compute_flux(G, np, trans, pmobf, gravcap_f, cpress, fflux);
|
||||||
|
|
||||||
if ((forces != NULL) && (forces->wells != NULL) &&
|
if ((forces != NULL) &&
|
||||||
(wpress != NULL) && (wflux != NULL)) {
|
(forces->wells != NULL) &&
|
||||||
|
(forces->wells->W != NULL) && (wpress != NULL) && (wflux != NULL)) {
|
||||||
|
|
||||||
compute_wflux(np, forces->wells, pmobc, cpress, wpress, wflux);
|
compute_wflux(np, forces->wells, pmobc, cpress, wpress, wflux);
|
||||||
}
|
}
|
||||||
|
|||||||
@@ -37,12 +37,20 @@ namespace Opm
|
|||||||
if (wells) {
|
if (wells) {
|
||||||
const int nw = wells->number_of_wells;
|
const int nw = wells->number_of_wells;
|
||||||
bhp_.resize(nw);
|
bhp_.resize(nw);
|
||||||
// Initialize bhp to be pressure in first perforation cell.
|
// Initialize bhp to be target pressure
|
||||||
|
// if bhp-controlled well, otherwise set
|
||||||
|
// to pressure in first perforation cell.
|
||||||
for (int w = 0; w < nw; ++w) {
|
for (int w = 0; w < nw; ++w) {
|
||||||
const int cell = wells->well_cells[wells->well_connpos[w]];
|
const WellControls* ctrl = wells->ctrls[w];
|
||||||
bhp_[w] = state.pressure()[cell];
|
if (ctrl->type[ctrl->current] == BHP) {
|
||||||
|
bhp_[w] = ctrl->target[ctrl->current];
|
||||||
|
} else {
|
||||||
|
const int cell = wells->well_cells[wells->well_connpos[w]];
|
||||||
|
bhp_[w] = state.pressure()[cell];
|
||||||
|
}
|
||||||
}
|
}
|
||||||
perfrates_.resize(wells->well_connpos[nw]);
|
perfrates_.resize(wells->well_connpos[nw], 0.0);
|
||||||
|
perfpress_.resize(wells->well_connpos[nw], -1e100);
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
|
|
||||||
@@ -54,9 +62,14 @@ namespace Opm
|
|||||||
std::vector<double>& perfRates() { return perfrates_; }
|
std::vector<double>& perfRates() { return perfrates_; }
|
||||||
const std::vector<double>& perfRates() const { return perfrates_; }
|
const std::vector<double>& perfRates() const { return perfrates_; }
|
||||||
|
|
||||||
|
/// One pressure per well connection.
|
||||||
|
std::vector<double>& perfPress() { return perfpress_; }
|
||||||
|
const std::vector<double>& perfPress() const { return perfpress_; }
|
||||||
|
|
||||||
private:
|
private:
|
||||||
std::vector<double> bhp_;
|
std::vector<double> bhp_;
|
||||||
std::vector<double> perfrates_;
|
std::vector<double> perfrates_;
|
||||||
|
std::vector<double> perfpress_;
|
||||||
};
|
};
|
||||||
|
|
||||||
} // namespace Opm
|
} // namespace Opm
|
||||||
|
|||||||
@@ -152,7 +152,7 @@ namespace Opm
|
|||||||
B_cell = 1.0/tm.A_[np*np*cell + 0];
|
B_cell = 1.0/tm.A_[np*np*cell + 0];
|
||||||
double src_flux = -tm.source_[cell];
|
double src_flux = -tm.source_[cell];
|
||||||
bool src_is_inflow = src_flux < 0.0;
|
bool src_is_inflow = src_flux < 0.0;
|
||||||
influx = src_is_inflow ? src_flux : 0.0;
|
influx = src_is_inflow ? B_cell* src_flux : 0.0;
|
||||||
outflux = !src_is_inflow ? src_flux : 0.0;
|
outflux = !src_is_inflow ? src_flux : 0.0;
|
||||||
comp_term = (tm.porevolume_[cell] - tm.porevolume0_[cell])/tm.porevolume0_[cell];
|
comp_term = (tm.porevolume_[cell] - tm.porevolume0_[cell])/tm.porevolume0_[cell];
|
||||||
dtpv = tm.dt_/tm.porevolume0_[cell];
|
dtpv = tm.dt_/tm.porevolume0_[cell];
|
||||||
|
|||||||
@@ -29,6 +29,7 @@
|
|||||||
#include <opm/core/utility/Units.hpp>
|
#include <opm/core/utility/Units.hpp>
|
||||||
#include <opm/core/fluid/IncompPropertiesInterface.hpp>
|
#include <opm/core/fluid/IncompPropertiesInterface.hpp>
|
||||||
#include <opm/core/fluid/BlackoilPropertiesInterface.hpp>
|
#include <opm/core/fluid/BlackoilPropertiesInterface.hpp>
|
||||||
|
#include <opm/core/fluid/blackoil/phaseUsageFromDeck.hpp>
|
||||||
#include <cmath>
|
#include <cmath>
|
||||||
|
|
||||||
namespace Opm
|
namespace Opm
|
||||||
@@ -512,15 +513,23 @@ namespace Opm
|
|||||||
State& state)
|
State& state)
|
||||||
{
|
{
|
||||||
const int num_phases = props.numPhases();
|
const int num_phases = props.numPhases();
|
||||||
|
const PhaseUsage pu = phaseUsageFromDeck(deck);
|
||||||
|
if (num_phases != pu.num_phases) {
|
||||||
|
THROW("initStateFromDeck(): user specified property object with " << num_phases << " phases, "
|
||||||
|
"found " << pu.num_phases << " phases in deck.");
|
||||||
|
}
|
||||||
state.init(grid, num_phases);
|
state.init(grid, num_phases);
|
||||||
if (deck.hasField("EQUIL")) {
|
if (deck.hasField("EQUIL")) {
|
||||||
if (num_phases != 2) {
|
if (num_phases != 2) {
|
||||||
THROW("initStateFromDeck(): EQUIL-based init currently handling only two-phase scenarios.");
|
THROW("initStateFromDeck(): EQUIL-based init currently handling only two-phase scenarios.");
|
||||||
}
|
}
|
||||||
|
if (pu.phase_used[BlackoilPhases::Vapour]) {
|
||||||
|
THROW("initStateFromDeck(): EQUIL-based init currently handling only oil-water scenario (no gas).");
|
||||||
|
}
|
||||||
// Set saturations depending on oil-water contact.
|
// Set saturations depending on oil-water contact.
|
||||||
const EQUIL& equil= deck.getEQUIL();
|
const EQUIL& equil= deck.getEQUIL();
|
||||||
if (equil.equil.size() != 1) {
|
if (equil.equil.size() != 1) {
|
||||||
THROW("No region support yet.");
|
THROW("initStateFromDeck(): No region support yet.");
|
||||||
}
|
}
|
||||||
const double woc = equil.equil[0].water_oil_contact_depth_;
|
const double woc = equil.equil[0].water_oil_contact_depth_;
|
||||||
initWaterOilContact(grid, props, woc, WaterBelow, state);
|
initWaterOilContact(grid, props, woc, WaterBelow, state);
|
||||||
@@ -528,37 +537,64 @@ namespace Opm
|
|||||||
const double datum_z = equil.equil[0].datum_depth_;
|
const double datum_z = equil.equil[0].datum_depth_;
|
||||||
const double datum_p = equil.equil[0].datum_depth_pressure_;
|
const double datum_p = equil.equil[0].datum_depth_pressure_;
|
||||||
initHydrostaticPressure(grid, props, woc, gravity, datum_z, datum_p, state);
|
initHydrostaticPressure(grid, props, woc, gravity, datum_z, datum_p, state);
|
||||||
} else if (deck.hasField("SWAT") && deck.hasField("PRESSURE")) {
|
} else if (deck.hasField("PRESSURE")) {
|
||||||
// Set saturations from SWAT, pressure from PRESSURE.
|
// Set saturations from SWAT/SGAS, pressure from PRESSURE.
|
||||||
std::vector<double>& s = state.saturation();
|
std::vector<double>& s = state.saturation();
|
||||||
std::vector<double>& p = state.pressure();
|
std::vector<double>& p = state.pressure();
|
||||||
const std::vector<double>& sw_deck = deck.getFloatingPointValue("SWAT");
|
|
||||||
const std::vector<double>& p_deck = deck.getFloatingPointValue("PRESSURE");
|
const std::vector<double>& p_deck = deck.getFloatingPointValue("PRESSURE");
|
||||||
const int num_cells = grid.number_of_cells;
|
const int num_cells = grid.number_of_cells;
|
||||||
if (num_phases == 2) {
|
if (num_phases == 2) {
|
||||||
for (int c = 0; c < num_cells; ++c) {
|
if (!pu.phase_used[BlackoilPhases::Aqua]) {
|
||||||
int c_deck = (grid.global_cell == NULL) ? c : grid.global_cell[c];
|
// oil-gas: we require SGAS
|
||||||
s[2*c] = sw_deck[c_deck];
|
if (!deck.hasField("SGAS")) {
|
||||||
s[2*c + 1] = 1.0 - s[2*c];
|
THROW("initStateFromDeck(): missing SGAS keyword in 2-phase init");
|
||||||
p[c] = p_deck[c_deck];
|
}
|
||||||
|
const std::vector<double>& sg_deck = deck.getFloatingPointValue("SGAS");
|
||||||
|
const int gpos = pu.phase_pos[BlackoilPhases::Vapour];
|
||||||
|
const int opos = pu.phase_pos[BlackoilPhases::Liquid];
|
||||||
|
for (int c = 0; c < num_cells; ++c) {
|
||||||
|
int c_deck = (grid.global_cell == NULL) ? c : grid.global_cell[c];
|
||||||
|
s[2*c + gpos] = sg_deck[c_deck];
|
||||||
|
s[2*c + opos] = 1.0 - sg_deck[c_deck];
|
||||||
|
p[c] = p_deck[c_deck];
|
||||||
|
}
|
||||||
|
} else {
|
||||||
|
// water-oil or water-gas: we require SWAT
|
||||||
|
if (!deck.hasField("SWAT")) {
|
||||||
|
THROW("initStateFromDeck(): missing SWAT keyword in 2-phase init");
|
||||||
|
}
|
||||||
|
const std::vector<double>& sw_deck = deck.getFloatingPointValue("SWAT");
|
||||||
|
const int wpos = pu.phase_pos[BlackoilPhases::Aqua];
|
||||||
|
const int nwpos = (wpos + 1) % 2;
|
||||||
|
for (int c = 0; c < num_cells; ++c) {
|
||||||
|
int c_deck = (grid.global_cell == NULL) ? c : grid.global_cell[c];
|
||||||
|
s[2*c + wpos] = sw_deck[c_deck];
|
||||||
|
s[2*c + nwpos] = 1.0 - sw_deck[c_deck];
|
||||||
|
p[c] = p_deck[c_deck];
|
||||||
|
}
|
||||||
}
|
}
|
||||||
} else if (num_phases == 3) {
|
} else if (num_phases == 3) {
|
||||||
if (!deck.hasField("SGAS")) {
|
const bool has_swat_sgas = deck.hasField("SWAT") && deck.hasField("SGAS");
|
||||||
THROW("initStateFromDeck(): missing SGAS keyword in 3-phase init (only SWAT and PRESSURE found).");
|
if (!has_swat_sgas) {
|
||||||
|
THROW("initStateFromDeck(): missing SGAS or SWAT keyword in 3-phase init.");
|
||||||
}
|
}
|
||||||
|
const int wpos = pu.phase_pos[BlackoilPhases::Aqua];
|
||||||
|
const int gpos = pu.phase_pos[BlackoilPhases::Vapour];
|
||||||
|
const int opos = pu.phase_pos[BlackoilPhases::Liquid];
|
||||||
|
const std::vector<double>& sw_deck = deck.getFloatingPointValue("SWAT");
|
||||||
const std::vector<double>& sg_deck = deck.getFloatingPointValue("SGAS");
|
const std::vector<double>& sg_deck = deck.getFloatingPointValue("SGAS");
|
||||||
for (int c = 0; c < num_cells; ++c) {
|
for (int c = 0; c < num_cells; ++c) {
|
||||||
int c_deck = (grid.global_cell == NULL) ? c : grid.global_cell[c];
|
int c_deck = (grid.global_cell == NULL) ? c : grid.global_cell[c];
|
||||||
s[3*c] = sw_deck[c_deck];
|
s[3*c + wpos] = sw_deck[c_deck];
|
||||||
s[3*c + 1] = 1.0 - (sw_deck[c_deck] + sg_deck[c_deck]);
|
s[3*c + opos] = 1.0 - (sw_deck[c_deck] + sg_deck[c_deck]);
|
||||||
s[3*c + 2] = sg_deck[c_deck];
|
s[3*c + gpos] = sg_deck[c_deck];
|
||||||
p[c] = p_deck[c_deck];
|
p[c] = p_deck[c_deck];
|
||||||
}
|
}
|
||||||
} else {
|
} else {
|
||||||
THROW("initStateFromDeck(): init with SWAT etc. only available with 2 or 3 phases.");
|
THROW("initStateFromDeck(): init with SWAT etc. only available with 2 or 3 phases.");
|
||||||
}
|
}
|
||||||
} else {
|
} else {
|
||||||
THROW("initStateFromDeck(): we must either have EQUIL, or both SWAT and PRESSURE.");
|
THROW("initStateFromDeck(): we must either have EQUIL, or PRESSURE and SWAT/SOIL/SGAS.");
|
||||||
}
|
}
|
||||||
|
|
||||||
// Finally, init face pressures.
|
// Finally, init face pressures.
|
||||||
|
|||||||
@@ -374,7 +374,10 @@ namespace Opm
|
|||||||
if (perf_rate > 0.0) {
|
if (perf_rate > 0.0) {
|
||||||
// perf_rate is a total inflow rate, we want a water rate.
|
// perf_rate is a total inflow rate, we want a water rate.
|
||||||
if (wells->type[w] != INJECTOR) {
|
if (wells->type[w] != INJECTOR) {
|
||||||
std::cout << "**** Warning: crossflow in well with index " << w << " ignored." << std::endl;
|
std::cout << "**** Warning: crossflow in well "
|
||||||
|
<< w << " perf " << perf - wells->well_connpos[w]
|
||||||
|
<< " ignored. Rate was "
|
||||||
|
<< perf_rate/Opm::unit::day << " m^3/day." << std::endl;
|
||||||
perf_rate = 0.0;
|
perf_rate = 0.0;
|
||||||
} else {
|
} else {
|
||||||
ASSERT(std::fabs(comp_frac[0] + comp_frac[1] - 1.0) < 1e-6);
|
ASSERT(std::fabs(comp_frac[0] + comp_frac[1] - 1.0) < 1e-6);
|
||||||
|
|||||||
@@ -21,7 +21,10 @@
|
|||||||
#include <opm/core/utility/miscUtilitiesBlackoil.hpp>
|
#include <opm/core/utility/miscUtilitiesBlackoil.hpp>
|
||||||
#include <opm/core/utility/Units.hpp>
|
#include <opm/core/utility/Units.hpp>
|
||||||
#include <opm/core/grid.h>
|
#include <opm/core/grid.h>
|
||||||
|
#include <opm/core/newwells.h>
|
||||||
#include <opm/core/fluid/BlackoilPropertiesInterface.hpp>
|
#include <opm/core/fluid/BlackoilPropertiesInterface.hpp>
|
||||||
|
#include <opm/core/simulator/BlackoilState.hpp>
|
||||||
|
#include <opm/core/simulator/WellState.hpp>
|
||||||
#include <opm/core/utility/ErrorMacros.hpp>
|
#include <opm/core/utility/ErrorMacros.hpp>
|
||||||
#include <algorithm>
|
#include <algorithm>
|
||||||
#include <functional>
|
#include <functional>
|
||||||
@@ -31,53 +34,74 @@
|
|||||||
namespace Opm
|
namespace Opm
|
||||||
{
|
{
|
||||||
|
|
||||||
/// @brief Computes injected and produced volumes of all phases.
|
/// @brief Computes injected and produced surface volumes of all phases.
|
||||||
/// Note 1: assumes that only the first phase is injected.
|
/// Note 1: assumes that only the first phase is injected.
|
||||||
/// Note 2: assumes that transport has been done with an
|
/// Note 2: assumes that transport has been done with an
|
||||||
/// implicit method, i.e. that the current state
|
/// implicit method, i.e. that the current state
|
||||||
/// gives the mobilities used for the preceding timestep.
|
/// gives the mobilities used for the preceding timestep.
|
||||||
/// @param[in] props fluid and rock properties.
|
/// Note 3: Gives surface volume values, not reservoir volumes
|
||||||
/// @param[in] p pressure (one value per cell)
|
/// (as the incompressible version of the function does).
|
||||||
/// @param[in] z surface-volume values (for all P phases)
|
/// Also, assumes that transport_src is given in surface volumes
|
||||||
/// @param[in] s saturation values (for all P phases)
|
/// for injector terms!
|
||||||
/// @param[in] src if < 0: total outflow, if > 0: first phase inflow.
|
/// @param[in] props fluid and rock properties.
|
||||||
/// @param[in] dt timestep used
|
/// @param[in] state state variables (pressure, sat, surfvol)
|
||||||
/// @param[out] injected must point to a valid array with P elements,
|
/// @param[in] transport_src if < 0: total resv outflow, if > 0: first phase surfv inflow
|
||||||
/// where P = s.size()/src.size().
|
/// @param[in] dt timestep used
|
||||||
/// @param[out] produced must also point to a valid array with P elements.
|
/// @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 BlackoilPropertiesInterface& props,
|
void computeInjectedProduced(const BlackoilPropertiesInterface& props,
|
||||||
const std::vector<double>& press,
|
const BlackoilState& state,
|
||||||
const std::vector<double>& z,
|
const std::vector<double>& transport_src,
|
||||||
const std::vector<double>& s,
|
|
||||||
const std::vector<double>& src,
|
|
||||||
const double dt,
|
const double dt,
|
||||||
double* injected,
|
double* injected,
|
||||||
double* produced)
|
double* produced)
|
||||||
{
|
{
|
||||||
const int num_cells = src.size();
|
const int num_cells = transport_src.size();
|
||||||
const int np = s.size()/src.size();
|
if (props.numCells() != num_cells) {
|
||||||
if (int(s.size()) != num_cells*np) {
|
THROW("Size of transport_src vector does not match number of cells in props.");
|
||||||
THROW("Sizes of s and src vectors do not match.");
|
|
||||||
}
|
}
|
||||||
|
const int np = props.numPhases();
|
||||||
|
if (int(state.saturation().size()) != num_cells*np) {
|
||||||
|
THROW("Sizes of state vectors do not match number of cells.");
|
||||||
|
}
|
||||||
|
const std::vector<double>& press = state.pressure();
|
||||||
|
const std::vector<double>& s = state.saturation();
|
||||||
|
const std::vector<double>& z = state.surfacevol();
|
||||||
std::fill(injected, injected + np, 0.0);
|
std::fill(injected, injected + np, 0.0);
|
||||||
std::fill(produced, produced + np, 0.0);
|
std::fill(produced, produced + np, 0.0);
|
||||||
std::vector<double> visc(np);
|
std::vector<double> visc(np);
|
||||||
std::vector<double> mob(np);
|
std::vector<double> mob(np);
|
||||||
|
std::vector<double> A(np*np);
|
||||||
|
std::vector<double> prod_resv_phase(np);
|
||||||
|
std::vector<double> prod_surfvol(np);
|
||||||
for (int c = 0; c < num_cells; ++c) {
|
for (int c = 0; c < num_cells; ++c) {
|
||||||
if (src[c] > 0.0) {
|
if (transport_src[c] > 0.0) {
|
||||||
injected[0] += src[c]*dt;
|
// Inflowing transport source is a surface volume flux
|
||||||
} else if (src[c] < 0.0) {
|
// for the first phase.
|
||||||
const double flux = -src[c]*dt;
|
injected[0] += transport_src[c]*dt;
|
||||||
|
} else if (transport_src[c] < 0.0) {
|
||||||
|
// Outflowing transport source is a total reservoir
|
||||||
|
// volume flux.
|
||||||
|
const double flux = -transport_src[c]*dt;
|
||||||
const double* sat = &s[np*c];
|
const double* sat = &s[np*c];
|
||||||
props.relperm(1, sat, &c, &mob[0], 0);
|
props.relperm(1, sat, &c, &mob[0], 0);
|
||||||
props.viscosity(1, &press[c], &z[np*c], &c, &visc[0], 0);
|
props.viscosity(1, &press[c], &z[np*c], &c, &visc[0], 0);
|
||||||
|
props.matrix(1, &press[c], &z[np*c], &c, &A[0], 0);
|
||||||
double totmob = 0.0;
|
double totmob = 0.0;
|
||||||
for (int p = 0; p < np; ++p) {
|
for (int p = 0; p < np; ++p) {
|
||||||
mob[p] /= visc[p];
|
mob[p] /= visc[p];
|
||||||
totmob += mob[p];
|
totmob += mob[p];
|
||||||
}
|
}
|
||||||
|
std::fill(prod_surfvol.begin(), prod_surfvol.end(), 0.0);
|
||||||
for (int p = 0; p < np; ++p) {
|
for (int p = 0; p < np; ++p) {
|
||||||
produced[p] += (mob[p]/totmob)*flux;
|
prod_resv_phase[p] = (mob[p]/totmob)*flux;
|
||||||
|
for (int q = 0; q < np; ++q) {
|
||||||
|
prod_surfvol[q] += prod_resv_phase[p]*A[q + np*p];
|
||||||
|
}
|
||||||
|
}
|
||||||
|
for (int p = 0; p < np; ++p) {
|
||||||
|
produced[p] += prod_surfvol[p];
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
@@ -251,4 +275,58 @@ namespace Opm
|
|||||||
}
|
}
|
||||||
}
|
}
|
||||||
|
|
||||||
|
|
||||||
|
/// Compute two-phase transport source terms from well terms.
|
||||||
|
/// Note: Unlike the incompressible version of this function,
|
||||||
|
/// this version computes surface volume injection rates,
|
||||||
|
/// production rates are still total reservoir volumes.
|
||||||
|
/// \param[in] props Fluid and rock properties.
|
||||||
|
/// \param[in] wells Wells data structure.
|
||||||
|
/// \param[in] well_state Well pressures and fluxes.
|
||||||
|
/// \param[out] transport_src The transport source terms. They are to be interpreted depending on sign:
|
||||||
|
/// (+) positive inflow of first (water) phase (surface volume),
|
||||||
|
/// (-) negative total outflow of both phases (reservoir volume).
|
||||||
|
void computeTransportSource(const BlackoilPropertiesInterface& props,
|
||||||
|
const Wells* wells,
|
||||||
|
const WellState& well_state,
|
||||||
|
std::vector<double>& transport_src)
|
||||||
|
{
|
||||||
|
int nc = props.numCells();
|
||||||
|
transport_src.clear();
|
||||||
|
transport_src.resize(nc, 0.0);
|
||||||
|
// Well contributions.
|
||||||
|
if (wells) {
|
||||||
|
const int nw = wells->number_of_wells;
|
||||||
|
const int np = wells->number_of_phases;
|
||||||
|
if (np != 2) {
|
||||||
|
THROW("computeTransportSource() requires a 2 phase case.");
|
||||||
|
}
|
||||||
|
std::vector<double> A(np*np);
|
||||||
|
for (int w = 0; w < nw; ++w) {
|
||||||
|
const double* comp_frac = wells->comp_frac + np*w;
|
||||||
|
for (int perf = wells->well_connpos[w]; perf < wells->well_connpos[w + 1]; ++perf) {
|
||||||
|
const int perf_cell = wells->well_cells[perf];
|
||||||
|
double perf_rate = well_state.perfRates()[perf];
|
||||||
|
if (perf_rate > 0.0) {
|
||||||
|
// perf_rate is a total inflow reservoir rate, we want a surface water rate.
|
||||||
|
if (wells->type[w] != INJECTOR) {
|
||||||
|
std::cout << "**** Warning: crossflow in well "
|
||||||
|
<< w << " perf " << perf - wells->well_connpos[w]
|
||||||
|
<< " ignored. Reservoir rate was "
|
||||||
|
<< perf_rate/Opm::unit::day << " m^3/day." << std::endl;
|
||||||
|
perf_rate = 0.0;
|
||||||
|
} else {
|
||||||
|
ASSERT(std::fabs(comp_frac[0] + comp_frac[1] - 1.0) < 1e-6);
|
||||||
|
perf_rate *= comp_frac[0]; // Water reservoir volume rate.
|
||||||
|
props.matrix(1, &well_state.perfPress()[perf], comp_frac, &perf_cell, &A[0], 0);
|
||||||
|
perf_rate *= A[0]; // Water surface volume rate.
|
||||||
|
}
|
||||||
|
}
|
||||||
|
transport_src[perf_cell] += perf_rate;
|
||||||
|
}
|
||||||
|
}
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
|
||||||
} // namespace Opm
|
} // namespace Opm
|
||||||
|
|||||||
@@ -22,36 +22,40 @@
|
|||||||
|
|
||||||
#include <vector>
|
#include <vector>
|
||||||
|
|
||||||
struct UnstructuredGrid;
|
struct Wells;
|
||||||
|
|
||||||
namespace Opm
|
namespace Opm
|
||||||
{
|
{
|
||||||
|
|
||||||
class BlackoilPropertiesInterface;
|
class BlackoilPropertiesInterface;
|
||||||
|
class BlackoilState;
|
||||||
|
class WellState;
|
||||||
|
|
||||||
/// @brief Computes injected and produced volumes of all phases.
|
|
||||||
|
/// @brief Computes injected and produced surface volumes of all phases.
|
||||||
/// Note 1: assumes that only the first phase is injected.
|
/// Note 1: assumes that only the first phase is injected.
|
||||||
/// Note 2: assumes that transport has been done with an
|
/// Note 2: assumes that transport has been done with an
|
||||||
/// implicit method, i.e. that the current state
|
/// implicit method, i.e. that the current state
|
||||||
/// gives the mobilities used for the preceding timestep.
|
/// gives the mobilities used for the preceding timestep.
|
||||||
/// @param[in] props fluid and rock properties.
|
/// Note 3: Gives surface volume values, not reservoir volumes
|
||||||
/// @param[in] p pressure (one value per cell)
|
/// (as the incompressible version of the function does).
|
||||||
/// @param[in] z surface-volume values (for all P phases)
|
/// Also, assumes that transport_src is given in surface volumes
|
||||||
/// @param[in] s saturation values (for all P phases)
|
/// for injector terms!
|
||||||
/// @param[in] src if < 0: total outflow, if > 0: first phase inflow.
|
/// @param[in] props fluid and rock properties.
|
||||||
/// @param[in] dt timestep used
|
/// @param[in] state state variables (pressure, sat, surfvol)
|
||||||
/// @param[out] injected must point to a valid array with P elements,
|
/// @param[in] transport_src if < 0: total resv outflow, if > 0: first phase surfv inflow
|
||||||
/// where P = s.size()/src.size().
|
/// @param[in] dt timestep used
|
||||||
/// @param[out] produced must also point to a valid array with P elements.
|
/// @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 BlackoilPropertiesInterface& props,
|
void computeInjectedProduced(const BlackoilPropertiesInterface& props,
|
||||||
const std::vector<double>& p,
|
const BlackoilState& state,
|
||||||
const std::vector<double>& z,
|
const std::vector<double>& transport_src,
|
||||||
const std::vector<double>& s,
|
|
||||||
const std::vector<double>& src,
|
|
||||||
const double dt,
|
const double dt,
|
||||||
double* injected,
|
double* injected,
|
||||||
double* produced);
|
double* produced);
|
||||||
|
|
||||||
|
|
||||||
/// @brief Computes total mobility for a set of saturation values.
|
/// @brief Computes total mobility for a set of saturation values.
|
||||||
/// @param[in] props rock and fluid properties
|
/// @param[in] props rock and fluid properties
|
||||||
/// @param[in] cells cells with which the saturation values are associated
|
/// @param[in] cells cells with which the saturation values are associated
|
||||||
@@ -66,6 +70,7 @@ namespace Opm
|
|||||||
const std::vector<double>& s,
|
const std::vector<double>& s,
|
||||||
std::vector<double>& totmob);
|
std::vector<double>& totmob);
|
||||||
|
|
||||||
|
|
||||||
/// @brief Computes total mobility and omega for a set of saturation values.
|
/// @brief Computes total mobility and omega for a set of saturation values.
|
||||||
/// @param[in] props rock and fluid properties
|
/// @param[in] props rock and fluid properties
|
||||||
/// @param[in] cells cells with which the saturation values are associated
|
/// @param[in] cells cells with which the saturation values are associated
|
||||||
@@ -131,6 +136,22 @@ namespace Opm
|
|||||||
const double* saturation,
|
const double* saturation,
|
||||||
double* surfacevol);
|
double* surfacevol);
|
||||||
|
|
||||||
|
|
||||||
|
/// Compute two-phase transport source terms from well terms.
|
||||||
|
/// Note: Unlike the incompressible version of this function,
|
||||||
|
/// this version computes surface volume injection rates,
|
||||||
|
/// production rates are still total reservoir volumes.
|
||||||
|
/// \param[in] props Fluid and rock properties.
|
||||||
|
/// \param[in] wells Wells data structure.
|
||||||
|
/// \param[in] well_state Well pressures and fluxes.
|
||||||
|
/// \param[out] transport_src The transport source terms. They are to be interpreted depending on sign:
|
||||||
|
/// (+) positive inflow of first (water) phase (surface volume),
|
||||||
|
/// (-) negative total outflow of both phases (reservoir volume).
|
||||||
|
void computeTransportSource(const BlackoilPropertiesInterface& props,
|
||||||
|
const Wells* wells,
|
||||||
|
const WellState& well_state,
|
||||||
|
std::vector<double>& transport_src);
|
||||||
|
|
||||||
} // namespace Opm
|
} // namespace Opm
|
||||||
|
|
||||||
#endif // OPM_MISCUTILITIESBLACKOIL_HEADER_INCLUDED
|
#endif // OPM_MISCUTILITIESBLACKOIL_HEADER_INCLUDED
|
||||||
|
|||||||
@@ -544,7 +544,7 @@ namespace Opm
|
|||||||
cf[pu.phase_pos[BlackoilPhases::Liquid]] = 1.0;
|
cf[pu.phase_pos[BlackoilPhases::Liquid]] = 1.0;
|
||||||
} else if (wci_line.injector_type_[0] == 'G') {
|
} else if (wci_line.injector_type_[0] == 'G') {
|
||||||
if (!pu.phase_used[BlackoilPhases::Vapour]) {
|
if (!pu.phase_used[BlackoilPhases::Vapour]) {
|
||||||
THROW("Water phase not used, yet found water-injecting well.");
|
THROW("Gas phase not used, yet found gas-injecting well.");
|
||||||
}
|
}
|
||||||
cf[pu.phase_pos[BlackoilPhases::Vapour]] = 1.0;
|
cf[pu.phase_pos[BlackoilPhases::Vapour]] = 1.0;
|
||||||
}
|
}
|
||||||
@@ -716,32 +716,31 @@ namespace Opm
|
|||||||
#endif
|
#endif
|
||||||
|
|
||||||
if (deck.hasField("WELOPEN")) {
|
if (deck.hasField("WELOPEN")) {
|
||||||
const WELOPEN& welopen = deck.getWELOPEN();
|
const WELOPEN& welopen = deck.getWELOPEN();
|
||||||
|
for (size_t i = 0; i < welopen.welopen.size(); ++i) {
|
||||||
for (size_t i = 0; i < welopen.welopen.size(); ++i) {
|
WelopenLine line = welopen.welopen[i];
|
||||||
WelopenLine line = welopen.welopen[i];
|
std::string wellname = line.well_;
|
||||||
std::string wellname = line.well_;
|
std::map<std::string, int>::const_iterator it = well_names_to_index.find(wellname);
|
||||||
std::map<std::string, int>::const_iterator it = well_names_to_index.find(wellname);
|
if (it == well_names_to_index.end()) {
|
||||||
if (it == well_names_to_index.end()) {
|
THROW("Trying to open/shut well with name: \"" << wellname<<"\" but it's not registered under WELSPECS.");
|
||||||
THROW("Trying to open/shut well with name: \"" << wellname<<"\" but it's not registered under WELSPECS.");
|
}
|
||||||
}
|
const int index = it->second;
|
||||||
int index = it->second;
|
if (line.openshutflag_ == "SHUT") {
|
||||||
if (line.openshutflag_ == "SHUT") {
|
int& cur_ctrl = w_->ctrls[index]->current;
|
||||||
// We currently don't care if the well is open or not.
|
if (cur_ctrl >= 0) {
|
||||||
/// \TODO Should this perhaps be allowed? I.e. should it be if(well_shut) { shutwell(); } else { /* do nothing*/ }?
|
cur_ctrl = ~cur_ctrl;
|
||||||
//ASSERT(w_->ctrls[index]->current < 0);
|
}
|
||||||
} else if (line.openshutflag_ == "OPEN") {
|
ASSERT(w_->ctrls[index]->current < 0);
|
||||||
//ASSERT(w_->ctrls[index]->current >= 0);
|
} else if (line.openshutflag_ == "OPEN") {
|
||||||
} else {
|
int& cur_ctrl = w_->ctrls[index]->current;
|
||||||
THROW("Unknown Open/close keyword: \"" << line.openshutflag_<< "\". Allowed values: OPEN, SHUT.");
|
if (cur_ctrl < 0) {
|
||||||
}
|
cur_ctrl = ~cur_ctrl;
|
||||||
|
}
|
||||||
// We revert back to it's original control.
|
ASSERT(w_->ctrls[index]->current >= 0);
|
||||||
// Note that this is OK as ~~ = id.
|
} else {
|
||||||
w_->ctrls[index]->current = ~w_->ctrls[index]->current;
|
THROW("Unknown Open/close keyword: \"" << line.openshutflag_<< "\". Allowed values: OPEN, SHUT.");
|
||||||
|
}
|
||||||
|
}
|
||||||
}
|
|
||||||
}
|
}
|
||||||
|
|
||||||
// Build the well_collection_ well group hierarchy.
|
// Build the well_collection_ well group hierarchy.
|
||||||
|
|||||||
Reference in New Issue
Block a user