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 pull request #430 from babrodtk/vfpprod

Browse Source 
Vertical flow performance
This commit is contained in:
Atgeirr Flø Rasmussen
2015-08-19 13:27:37 +02:00
parent dbc9aedc5a 407512382a
commit 3c1c3a30ab
20 changed files with 8160 additions and 66 deletions
Download Patch File Download Diff File Expand all files Collapse all files

10
CMakeLists_files.cmake
View File

@@ -42,6 +42,9 @@ list (APPEND MAIN_SOURCE_FILES
opm/autodiff/BlackoilModelParameters.cpp
opm/autodiff/WellDensitySegmented.cpp
opm/autodiff/LinearisedBlackoilResidual.cpp
opm/autodiff/VFPProperties.cpp
opm/autodiff/VFPProdProperties.cpp
opm/autodiff/VFPInjProperties.cpp
)
# originally generated with the command:
@@ -56,10 +59,13 @@ list (APPEND TEST_SOURCE_FILES
tests/test_scalar_mult.cpp
tests/test_transmissibilitymultipliers.cpp
tests/test_welldensitysegmented.cpp
tests/test_vfpproperties.cpp
)
list (APPEND TEST_DATA_FILES
tests/fluid.data
tests/VFPPROD1
tests/VFPPROD2
)
# Note, these two files are not included in the repo.
@@ -138,5 +144,9 @@ list (APPEND PUBLIC_HEADER_FILES
opm/autodiff/WellStateFullyImplicitBlackoil.hpp
opm/autodiff/WellStateFullyImplicitBlackoilSolvent.hpp
opm/autodiff/SimulatorFullyImplicitBlackoilOutput.hpp
opm/autodiff/VFPProperties.hpp
opm/autodiff/VFPHelpers.hpp
opm/autodiff/VFPProdProperties.hpp
opm/autodiff/VFPInjProperties.hpp
)

11
opm/autodiff/AutoDiffHelpers.hpp
View File

@@ -263,6 +263,12 @@ namespace {
constructSupersetSparseMatrix(const int full_size, const IntVec& indices)
{
const int subset_size = indices.size();
if (subset_size == 0) {
typename AutoDiffBlock<Scalar>::M mat(full_size, 0);
return mat;
}
typename AutoDiffBlock<Scalar>::M mat(full_size, subset_size);
mat.reserve(Eigen::VectorXi::Constant(subset_size, 1));
for (int i = 0; i < subset_size; ++i) {
@@ -354,7 +360,7 @@ spdiag(const AutoDiffBlock<double>::V& d)
public:
typedef AutoDiffBlock<Scalar> ADB;
enum CriterionForLeftElement { GreaterEqualZero, GreaterZero, Zero, NotEqualZero, LessZero, LessEqualZero };
enum CriterionForLeftElement { GreaterEqualZero, GreaterZero, Zero, NotEqualZero, LessZero, LessEqualZero, NotNaN };
Selector(const typename ADB::V& selection_basis,
CriterionForLeftElement crit = GreaterEqualZero)
@@ -385,6 +391,9 @@ spdiag(const AutoDiffBlock<double>::V& d)
case LessEqualZero:
chooseleft = selection_basis[i] <= 0.0;
break;
case NotNaN:
chooseleft = !isnan(selection_basis[i]);
break;
default:
OPM_THROW(std::logic_error, "No such criterion: " << crit);
}

7
opm/autodiff/BlackoilModelBase.hpp
View File

@@ -31,6 +31,7 @@
#include <opm/autodiff/LinearisedBlackoilResidual.hpp>
#include <opm/autodiff/NewtonIterationBlackoilInterface.hpp>
#include <opm/autodiff/BlackoilModelEnums.hpp>
#include <opm/autodiff/VFPProperties.hpp>
#include <opm/parser/eclipse/EclipseState/Grid/NNC.hpp>
#include <array>
@@ -43,6 +44,7 @@ namespace Opm {
class DerivedGeology;
class RockCompressibility;
class NewtonIterationBlackoilInterface;
class VFPProperties;
/// Struct for containing iteration variables.
@@ -117,6 +119,7 @@ namespace Opm {
/// \param[in] geo rock properties
/// \param[in] rock_comp_props if non-null, rock compressibility properties
/// \param[in] wells well structure
/// \param[in] vfp_properties Vertical flow performance tables
/// \param[in] linsolver linear solver
/// \param[in] eclState eclipse state
/// \param[in] has_disgas turn on dissolved gas
@@ -236,6 +239,7 @@ namespace Opm {
const DerivedGeology& geo_;
const RockCompressibility* rock_comp_props_;
const Wells* wells_;
VFPProperties vfp_properties_;
const NewtonIterationBlackoilInterface& linsolver_;
// For each canonical phase -> true if active
const std::vector<bool> active_;
@@ -256,6 +260,7 @@ namespace Opm {
V isRs_;
V isRv_;
V isSg_;
V well_perforation_densities_; //Density of each well perforation
V well_perforation_pressure_diffs_; // Diff to bhp for each well perforation.
LinearisedBlackoilResidual residual_;
@@ -369,6 +374,8 @@ namespace Opm {
bool getWellConvergence(const int iteration);
bool isVFPActive() const;
std::vector<ADB>
computePressures(const ADB& po,
const ADB& sw,

340
opm/autodiff/BlackoilModelBase_impl.hpp
View File

@@ -32,6 +32,9 @@
#include <opm/autodiff/BlackoilPropsAdInterface.hpp>
#include <opm/autodiff/GeoProps.hpp>
#include <opm/autodiff/WellDensitySegmented.hpp>
#include <opm/autodiff/VFPProperties.hpp>
#include <opm/autodiff/VFPProdProperties.hpp>
#include <opm/autodiff/VFPInjProperties.hpp>
#include <opm/core/grid.h>
#include <opm/core/linalg/LinearSolverInterface.hpp>
@@ -49,6 +52,7 @@
#include <iostream>
#include <iomanip>
#include <limits>
#include <vector>
//#include <fstream>
// A debugging utility.
@@ -132,7 +136,6 @@ namespace detail {
return act2can;
}
} // namespace detail
@@ -154,6 +157,7 @@ namespace detail {
, geo_ (geo)
, rock_comp_props_(rock_comp_props)
, wells_ (wells)
, vfp_properties_(eclState->getVFPInjTables(), eclState->getVFPProdTables())
, linsolver_ (linsolver)
, active_(detail::activePhases(fluid.phaseUsage()))
, canph_ (detail::active2Canonical(fluid.phaseUsage()))
@@ -652,7 +656,19 @@ namespace detail {
}
}
namespace detail {
double getGravity(const double* g, const int dim) {
double grav = 0.0;
if (g) {
// Guard against gravity in anything but last dimension.
for (int dd = 0; dd < dim - 1; ++dd) {
assert(g[dd] == 0.0);
}
grav = g[dim - 1];
}
return grav;
}
}
@@ -726,22 +742,21 @@ namespace detail {
// Surface density.
std::vector<double> surf_dens(fluid_.surfaceDensity(), fluid_.surfaceDensity() + pu.num_phases);
// Gravity
double grav = 0.0;
const double* g = geo_.gravity();
const int dim = dimensions(grid_);
if (g) {
// Guard against gravity in anything but last dimension.
for (int dd = 0; dd < dim - 1; ++dd) {
assert(g[dd] == 0.0);
}
grav = g[dim - 1];
}
double grav = detail::getGravity(geo_.gravity(), dimensions(grid_));
// 2. Compute pressure deltas, and store the results.
std::vector<double> cdp = WellDensitySegmented
::computeConnectionPressureDelta(wells(), xw, fluid_.phaseUsage(),
b_perf, rsmax_perf, rvmax_perf, perf_depth,
surf_dens, grav);
// 2. Compute densities
std::vector<double> cd =
WellDensitySegmented::computeConnectionDensities(
wells(), xw, fluid_.phaseUsage(),
b_perf, rsmax_perf, rvmax_perf, surf_dens);
// 3. Compute pressure deltas
std::vector<double> cdp =
WellDensitySegmented::computeConnectionPressureDelta(
wells(), perf_depth, cd, grav);
// 4. Store the results
well_perforation_densities_ = Eigen::Map<const V>(cd.data(), nperf);
well_perforation_pressure_diffs_ = Eigen::Map<const V>(cdp.data(), nperf);
}
@@ -758,6 +773,16 @@ namespace detail {
{
using namespace Opm::AutoDiffGrid;
// If we have VFP tables, we need the well connection
// pressures for the "simple" hydrostatic correction
// between well depth and vfp table depth.
if (isVFPActive()) {
SolutionState state = asImpl().variableState(reservoir_state, well_state);
SolutionState state0 = state;
asImpl().makeConstantState(state0);
computeWellConnectionPressures(state0, well_state);
}
// Possibly switch well controls and updating well state to
// get reasonable initial conditions for the wells
updateWellControls(well_state);
@@ -816,7 +841,7 @@ namespace detail {
asImpl().updatePerfPhaseRatesAndPressures(cq_s, state, well_state);
asImpl().addWellFluxEq(cq_s, state);
asImpl().addWellContributionToMassBalanceEq(cq_s, state, well_state);
addWellControlEq(state, well_state, aliveWells);
addWellControlEq(state, well_state, aliveWells);
}
@@ -901,7 +926,6 @@ namespace detail {
template <class Grid, class Implementation>
void
BlackoilModelBase<Grid, Implementation>::computeWellFlux(const SolutionState& state,
@@ -1094,6 +1118,7 @@ namespace detail {
}
bool constraintBroken(const std::vector<double>& bhp,
const std::vector<double>& thp,
const std::vector<double>& well_phase_flow_rate,
const int well,
const int num_phases,
@@ -1115,6 +1140,10 @@ namespace detail {
broken = bhp[well] > target;
break;
case THP:
broken = thp[well] > target;
break;
case RESERVOIR_RATE: // Intentional fall-through
case SURFACE_RATE:
broken = rateToCompare(well_phase_flow_rate,
@@ -1131,6 +1160,10 @@ namespace detail {
broken = bhp[well] < target;
break;
case THP:
broken = thp[well] < target;
break;
case RESERVOIR_RATE: // Intentional fall-through
case SURFACE_RATE:
// Note that the rates compared below are negative,
@@ -1152,7 +1185,51 @@ namespace detail {
} // namespace detail
namespace detail {
double computeHydrostaticCorrection(const Wells& wells, const int w, const double vfp_ref_depth,
const ADB::V& well_perforation_densities, const double gravity) {
//For the initial iteration, we have no perforation densities.
const double well_ref_depth = wells.depth_ref[w];
const double dh = vfp_ref_depth - well_ref_depth;
const int perf = wells.well_connpos[w];
const double rho = well_perforation_densities[perf];
const double dp = rho*gravity*dh;
return dp;
}
} //Namespace
template <class Grid, class Implementation>
bool BlackoilModelBase<Grid, Implementation>::isVFPActive() const
{
if( ! wellsActive() ) {
return false;
}
if ( vfp_properties_.getProd()->empty() && vfp_properties_.getInj()->empty() ) {
return false;
}
const int nw = wells().number_of_wells;
//Loop over all wells
for (int w = 0; w < nw; ++w) {
const WellControls* wc = wells().ctrls[w];
const int nwc = well_controls_get_num(wc);
//Loop over all controls
for (int c=0; c < nwc; ++c) {
const WellControlType ctrl_type = well_controls_iget_type(wc, c);
if (ctrl_type == THP) {
return true;
}
}
}
return false;
}
template <class Grid, class Implementation>
@@ -1160,11 +1237,12 @@ namespace detail {
{
if( ! wellsActive() ) return ;
std::string modestring[3] = { "BHP", "RESERVOIR_RATE", "SURFACE_RATE" };
std::string modestring[4] = { "BHP", "THP", "RESERVOIR_RATE", "SURFACE_RATE" };
// Find, for each well, if any constraints are broken. If so,
// switch control to first broken constraint.
const int np = wells().number_of_phases;
const int nw = wells().number_of_wells;
const Opm::PhaseUsage& pu = fluid_.phaseUsage();
for (int w = 0; w < nw; ++w) {
const WellControls* wc = wells().ctrls[w];
// The current control in the well state overrides
@@ -1183,7 +1261,9 @@ namespace detail {
// inequality constraint, and therefore skipped.
continue;
}
if (detail::constraintBroken(xw.bhp(), xw.wellRates(), w, np, wells().type[w], wc, ctrl_index)) {
if (detail::constraintBroken(
xw.bhp(), xw.thp(), xw.wellRates(),
w, np, wells().type[w], wc, ctrl_index)) {
// ctrl_index will be the index of the broken constraint after the loop.
break;
}
@@ -1200,8 +1280,11 @@ namespace detail {
current = xw.currentControls()[w];
}
//Get gravity for THP hydrostatic corrrection
const double gravity = detail::getGravity(geo_.gravity(), UgGridHelpers::dimensions(grid_));
// Updating well state and primary variables.
// Target values are used as initial conditions for BHP and SURFACE_RATE
// Target values are used as initial conditions for BHP, THP, and SURFACE_RATE
const double target = well_controls_iget_target(wc, current);
const double* distr = well_controls_iget_distr(wc, current);
switch (well_controls_iget_type(wc, current)) {
@@ -1209,6 +1292,52 @@ namespace detail {
xw.bhp()[w] = target;
break;
case THP: {
double aqua = 0.0;
double liquid = 0.0;
double vapour = 0.0;
if (active_[ Water ]) {
aqua = xw.wellRates()[w*np + pu.phase_pos[ Water ] ];
}
if (active_[ Oil ]) {
liquid = xw.wellRates()[w*np + pu.phase_pos[ Oil ] ];
}
if (active_[ Gas ]) {
vapour = xw.wellRates()[w*np + pu.phase_pos[ Gas ] ];
}
const int vfp = well_controls_iget_vfp(wc, current);
const double& thp = well_controls_iget_target(wc, current);
const double& alq = well_controls_iget_alq(wc, current);
//Set *BHP* target by calculating bhp from THP
const WellType& well_type = wells().type[w];
//Gather variables for hydrostatic reconstruction
if (well_type == INJECTOR) {
double vfp_ref_depth = vfp_properties_.getInj()->getTable(vfp)->getDatumDepth();
double dp = detail::computeHydrostaticCorrection(
wells(), w, vfp_ref_depth,
well_perforation_densities_, gravity);
xw.bhp()[w] = vfp_properties_.getInj()->bhp(vfp, aqua, liquid, vapour, thp) - dp;
}
else if (well_type == PRODUCER) {
double vfp_ref_depth = vfp_properties_.getProd()->getTable(vfp)->getDatumDepth();
double dp = detail::computeHydrostaticCorrection(
wells(), w, vfp_ref_depth,
well_perforation_densities_, gravity);
xw.bhp()[w] = vfp_properties_.getProd()->bhp(vfp, aqua, liquid, vapour, thp, alq) - dp;
}
else {
OPM_THROW(std::logic_error, "Expected PRODUCER or INJECTOR type of well");
}
break;
}
case RESERVOIR_RATE:
// No direct change to any observable quantity at
// surface condition. In this case, use existing
@@ -1328,11 +1457,47 @@ namespace detail {
const int np = wells().number_of_phases;
const int nw = wells().number_of_wells;
V bhp_targets = V::Zero(nw);
V rate_targets = V::Zero(nw);
M rate_distr(nw, np*nw);
ADB aqua = ADB::constant(ADB::V::Zero(nw));
ADB liquid = ADB::constant(ADB::V::Zero(nw));
ADB vapour = ADB::constant(ADB::V::Zero(nw));
if (active_[Water]) {
aqua += subset(state.qs, Span(nw, 1, BlackoilPhases::Aqua*nw));
}
if (active_[Oil]) {
liquid += subset(state.qs, Span(nw, 1, BlackoilPhases::Liquid*nw));
}
if (active_[Gas]) {
vapour += subset(state.qs, Span(nw, 1, BlackoilPhases::Vapour*nw));
}
//THP calculation variables
std::vector<int> inj_table_id(nw, -1);
std::vector<int> prod_table_id(nw, -1);
ADB::V thp_inj_target_v = ADB::V::Zero(nw);
ADB::V thp_prod_target_v = ADB::V::Zero(nw);
ADB::V alq_v = ADB::V::Zero(nw);
//Hydrostatic correction variables
ADB::V rho_v = ADB::V::Zero(nw);
ADB::V vfp_ref_depth_v = ADB::V::Zero(nw);
//Target vars
ADB::V bhp_targets = ADB::V::Zero(nw);
ADB::V rate_targets = ADB::V::Zero(nw);
ADB::M rate_distr(nw, np*nw);
//Selection variables
std::vector<int> bhp_elems;
std::vector<int> thp_inj_elems;
std::vector<int> thp_prod_elems;
std::vector<int> rate_elems;
//Run through all wells to calculate BHP/RATE targets
//and gather info about current control
for (int w = 0; w < nw; ++w) {
const WellControls* wc = wells().ctrls[w];
auto wc = wells().ctrls[w];
// The current control in the well state overrides
// the current control set in the Wells struct, which
// is instead treated as a default.
@@ -1341,7 +1506,43 @@ namespace detail {
switch (well_controls_iget_type(wc, current)) {
case BHP:
{
bhp_targets (w) = well_controls_iget_target(wc, current);
bhp_elems.push_back(w);
bhp_targets(w) = well_controls_iget_target(wc, current);
rate_targets(w) = -1e100;
}
break;
case THP:
{
const int perf = wells().well_connpos[w];
rho_v[w] = well_perforation_densities_[perf];
const int table_id = well_controls_iget_vfp(wc, current);
const double target = well_controls_iget_target(wc, current);
const WellType& well_type = wells().type[w];
if (well_type == INJECTOR) {
inj_table_id[w] = table_id;
thp_inj_target_v[w] = target;
alq_v[w] = -1e100;
vfp_ref_depth_v[w] = vfp_properties_.getInj()->getTable(table_id)->getDatumDepth();
thp_inj_elems.push_back(w);
}
else if (well_type == PRODUCER) {
prod_table_id[w] = table_id;
thp_prod_target_v[w] = target;
alq_v[w] = well_controls_iget_alq(wc, current);
vfp_ref_depth_v[w] = vfp_properties_.getProd()->getTable(table_id)->getDatumDepth();
thp_prod_elems.push_back(w);
}
else {
OPM_THROW(std::logic_error, "Expected INJECTOR or PRODUCER type well");
}
bhp_targets(w) = -1e100;
rate_targets(w) = -1e100;
}
break;
@@ -1349,6 +1550,7 @@ namespace detail {
case RESERVOIR_RATE: // Intentional fall-through
case SURFACE_RATE:
{
rate_elems.push_back(w);
// RESERVOIR and SURFACE rates look the same, from a
// high-level point of view, in the system of
// simultaneous linear equations.
@@ -1360,17 +1562,43 @@ namespace detail {
rate_distr.insert(w, p*nw + w) = distr[p];
}
bhp_targets (w) = -1.0e100;
bhp_targets(w) = -1.0e100;
rate_targets(w) = well_controls_iget_target(wc, current);
}
break;
}
}
//Calculate BHP target from THP
const ADB thp_inj_target = ADB::constant(thp_inj_target_v);
const ADB thp_prod_target = ADB::constant(thp_prod_target_v);
const ADB alq = ADB::constant(alq_v);
const ADB bhp_from_thp_inj = vfp_properties_.getInj()->bhp(inj_table_id, aqua, liquid, vapour, thp_inj_target);
const ADB bhp_from_thp_prod = vfp_properties_.getProd()->bhp(prod_table_id, aqua, liquid, vapour, thp_prod_target, alq);
//Perform hydrostatic correction to computed targets
//FIXME: Use computeHydrostaticCorrection
const ADB well_ref_depth = ADB::constant(ADB::V::Map(wells().depth_ref, nw));
const ADB vfp_ref_depth = ADB::constant(vfp_ref_depth_v);
const ADB dh = vfp_ref_depth - well_ref_depth;
const ADB rho = ADB::constant(rho_v);
double gravity = detail::getGravity(geo_.gravity(), UgGridHelpers::dimensions(grid_));
const ADB dp = rho*gravity*dh;
const ADB dp_inj = superset(subset(dp, thp_inj_elems), thp_inj_elems, nw);
const ADB dp_prod = superset(subset(dp, thp_prod_elems), thp_prod_elems, nw);
//Calculate residuals
const ADB thp_inj_residual = state.bhp - bhp_from_thp_inj + dp_inj;
const ADB thp_prod_residual = state.bhp - bhp_from_thp_prod + dp_prod;
const ADB bhp_residual = state.bhp - bhp_targets;
const ADB rate_residual = rate_distr * state.qs - rate_targets;
// Choose bhp residual for positive bhp targets.
Selector<double> bhp_selector(bhp_targets);
residual_.well_eq = bhp_selector.select(bhp_residual, rate_residual);
//Select the right residual for each well
residual_.well_eq = superset(subset(bhp_residual, bhp_elems), bhp_elems, nw) +
superset(subset(thp_inj_residual, thp_inj_elems), thp_inj_elems, nw) +
superset(subset(thp_prod_residual, thp_prod_elems), thp_prod_elems, nw) +
superset(subset(rate_residual, rate_elems), rate_elems, nw);
// For wells that are dead (not flowing), and therefore not communicating
// with the reservoir, we set the equation to be equal to the well's total
// flow. This will be a solution only if the target rate is also zero.
@@ -1716,6 +1944,51 @@ namespace detail {
const V dbhp_limited = sign(dbhp) * dbhp.abs().min(bhp_old.abs()*dpmaxrel);
const V bhp = bhp_old - dbhp_limited;
std::copy(&bhp[0], &bhp[0] + bhp.size(), well_state.bhp().begin());
// Thp update
const Opm::PhaseUsage& pu = fluid_.phaseUsage();
//Loop over all wells
for (int w=0; w<nw; ++w) {
const WellControls* wc = wells().ctrls[w];
const int nwc = well_controls_get_num(wc);
//Loop over all controls until we find a THP control
//that specifies what we need...
//Will only update THP for wells with THP control
for (int ctrl_index=0; ctrl_index < nwc; ++ctrl_index) {
if (well_controls_iget_type(wc, ctrl_index) == THP) {
double aqua = 0.0;
double liquid = 0.0;
double vapour = 0.0;
if (active_[ Water ]) {
aqua = wr[w*np + pu.phase_pos[ Water ] ];
}
if (active_[ Oil ]) {
liquid = wr[w*np + pu.phase_pos[ Oil ] ];
}
if (active_[ Gas ]) {
vapour = wr[w*np + pu.phase_pos[ Gas ] ];
}
double alq = well_controls_iget_alq(wc, ctrl_index);
int table_id = well_controls_iget_vfp(wc, ctrl_index);
const WellType& well_type = wells().type[w];
if (well_type == INJECTOR) {
well_state.thp()[w] = vfp_properties_.getInj()->thp(table_id, aqua, liquid, vapour, bhp[w]);
}
else if (well_type == PRODUCER) {
well_state.thp()[w] = vfp_properties_.getProd()->thp(table_id, aqua, liquid, vapour, bhp[w], alq);
}
else {
OPM_THROW(std::logic_error, "Expected INJECTOR or PRODUCER well");
}
//Assume only one THP control specified for each well
break;
}
}
}
}
}
@@ -2065,6 +2338,9 @@ namespace detail {
converged_Well = converged_Well && (residualWell < Opm::unit::barsa);
const bool converged = converged_MB && converged_CNV && converged_Well;
// Residual in Pascal can have high values and still be ok.
const double maxWellResidualAllowed = 1000.0 * maxResidualAllowed();
// if one of the residuals is NaN, throw exception, so that the solver can be restarted
if ( std::isnan(mass_balance_residual[Water]) || mass_balance_residual[Water] > maxResidualAllowed() ||
std::isnan(mass_balance_residual[Oil]) || mass_balance_residual[Oil] > maxResidualAllowed() ||
@@ -2075,7 +2351,7 @@ namespace detail {
std::isnan(well_flux_residual[Water]) || well_flux_residual[Water] > maxResidualAllowed() ||
std::isnan(well_flux_residual[Oil]) || well_flux_residual[Oil] > maxResidualAllowed() ||
std::isnan(well_flux_residual[Gas]) || well_flux_residual[Gas] > maxResidualAllowed() ||
std::isnan(residualWell) || residualWell > maxResidualAllowed() )
std::isnan(residualWell) || residualWell > maxWellResidualAllowed )
{
OPM_THROW(Opm::NumericalProblem,"One of the residuals is NaN or too large!");
}

2
opm/autodiff/BlackoilPropsAdFromDeck.hpp
View File

@@ -347,7 +347,7 @@ namespace Opm
// Densities, one std::array per PVT region.
std::vector<std::array<double, BlackoilPhases::MaxNumPhases> > densities_;
// VAPPARS
double vap1_;
double vap2_;

5
opm/autodiff/SimulatorBase_impl.hpp
View File

@@ -437,8 +437,12 @@ namespace Opm
well_controls_clear(ctrl);
well_controls_assert_number_of_phases(ctrl, int(np));
static const double invalid_alq = -std::numeric_limits<double>::max();
static const int invalid_vfp = -std::numeric_limits<int>::max();
const int ok_resv =
well_controls_add_new(RESERVOIR_RATE, target,
invalid_alq, invalid_vfp,
& distr[0], ctrl);
// For WCONHIST/RESV the BHP limit is set to 1 atm.
@@ -446,6 +450,7 @@ namespace Opm
// the WELTARG keyword
const int ok_bhp =
well_controls_add_new(BHP, unit::convert::from(1.0, unit::atm),
invalid_alq, invalid_vfp,
NULL, ctrl);
if (ok_resv != 0 && ok_bhp != 0) {

983
opm/autodiff/VFPHelpers.hpp Normal file
View File

@@ -0,0 +1,983 @@
/*
Copyright 2015 SINTEF ICT, Applied Mathematics.
This file is part of the Open Porous Media project (OPM).
OPM is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OPM is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef OPM_AUTODIFF_VFPHELPERS_HPP_
#define OPM_AUTODIFF_VFPHELPERS_HPP_
#include <opm/parser/eclipse/EclipseState/Tables/VFPProdTable.hpp>
#include <opm/parser/eclipse/EclipseState/Tables/VFPInjTable.hpp>
#include <opm/autodiff/AutoDiffHelpers.hpp>
/**
* This file contains a set of helper functions used by VFPProd / VFPInj.
*/
namespace Opm {
namespace detail {
typedef AutoDiffBlock<double> ADB;
/**
* Returns zero if input value is NaN
*/
inline double zeroIfNan(const double& value) {
return (std::isnan(value)) ? 0.0 : value;
}
/**
* Returns zero for every entry in the ADB which is NaN
*/
inline ADB zeroIfNan(const ADB& values) {
Selector<ADB::V::Scalar> not_nan_selector(values.value(), Selector<ADB::V::Scalar>::NotNaN);
const ADB::V z = ADB::V::Zero(values.size());
const ADB zero = ADB::constant(z, values.blockPattern());
ADB retval = not_nan_selector.select(values, zero);
return retval;
}
/**
* Computes the flo parameter according to the flo_type_
* for production tables
* @return Production rate of oil, gas or liquid.
*/
template <typename T>
static T getFlo(const T& aqua, const T& liquid, const T& vapour,
const VFPProdTable::FLO_TYPE& type) {
switch (type) {
case VFPProdTable::FLO_OIL:
//Oil = liquid phase
return liquid;
case VFPProdTable::FLO_LIQ:
//Liquid = aqua + liquid phases
return aqua + liquid;
case VFPProdTable::FLO_GAS:
//Gas = vapor phase
return vapour;
case VFPProdTable::FLO_INVALID: //Intentional fall-through
default:
OPM_THROW(std::logic_error, "Invalid FLO_TYPE: '" << type << "'");
}
}
/**
* Computes the flo parameter according to the flo_type_
* for injection tables
* @return Production rate of oil, gas or liquid.
*/
template <typename T>
static T getFlo(const T& aqua, const T& liquid, const T& vapour,
const VFPInjTable::FLO_TYPE& type) {
switch (type) {
case VFPInjTable::FLO_OIL:
//Oil = liquid phase
return liquid;
case VFPInjTable::FLO_WAT:
//Liquid = aqua phase
return aqua;
case VFPInjTable::FLO_GAS:
//Gas = vapor phase
return vapour;
case VFPInjTable::FLO_INVALID: //Intentional fall-through
default:
OPM_THROW(std::logic_error, "Invalid FLO_TYPE: '" << type << "'");
}
}
/**
* Computes the wfr parameter according to the wfr_type_
* @return Production rate of oil, gas or liquid.
*/
template <typename T>
static T getWFR(const T& aqua, const T& liquid, const T& vapour,
const VFPProdTable::WFR_TYPE& type) {
switch(type) {
case VFPProdTable::WFR_WOR: {
//Water-oil ratio = water / oil
T wor = aqua / liquid;
return zeroIfNan(wor);
}
case VFPProdTable::WFR_WCT:
//Water cut = water / (water + oil)
return zeroIfNan(aqua / (aqua + liquid));
case VFPProdTable::WFR_WGR:
//Water-gas ratio = water / gas
return zeroIfNan(aqua / vapour);
case VFPProdTable::WFR_INVALID: //Intentional fall-through
default:
OPM_THROW(std::logic_error, "Invalid WFR_TYPE: '" << type << "'");
}
}
/**
* Computes the gfr parameter according to the gfr_type_
* @return Production rate of oil, gas or liquid.
*/
template <typename T>
static T getGFR(const T& aqua, const T& liquid, const T& vapour,
const VFPProdTable::GFR_TYPE& type) {
switch(type) {
case VFPProdTable::GFR_GOR:
// Gas-oil ratio = gas / oil
return zeroIfNan(vapour / liquid);
case VFPProdTable::GFR_GLR:
// Gas-liquid ratio = gas / (oil + water)
return zeroIfNan(vapour / (liquid + aqua));
case VFPProdTable::GFR_OGR:
// Oil-gas ratio = oil / gas
return zeroIfNan(liquid / vapour);
case VFPProdTable::GFR_INVALID: //Intentional fall-through
default:
OPM_THROW(std::logic_error, "Invalid GFR_TYPE: '" << type << "'");
}
}
/**
* Helper struct for linear interpolation
*/
struct InterpData {
InterpData() : ind_{0, 0}, inv_dist_(0.0), factor_(0.0) {}
int ind_[2]; //[First element greater than or equal to value, Last element smaller than or equal to value]
double inv_dist_; // 1 / distance between the two end points of the segment. Used to calculate derivatives and uses 1.0 / 0.0 = 0.0 as a convention
double factor_; // Interpolation factor
};
/**
* Helper function to find indices etc. for linear interpolation and extrapolation
* @param value Value to find in values
* @param values Sorted list of values to search for value in.
* @return Data required to find the interpolated value
*/
inline InterpData findInterpData(const double& value, const std::vector<double>& values) {
InterpData retval;
const int nvalues = values.size();
//If we only have one value in our vector, return that
if (values.size() == 1) {
retval.ind_[0] = 0;
retval.ind_[1] = 0;
retval.inv_dist_ = 0.0;
retval.factor_ = 0.0;
}
// Else search in the vector
else {
//If value is less than all values, use first interval
if (value < values.front()) {
retval.ind_[0] = 0;
retval.ind_[1] = 1;
}
//If value is greater than all values, use last interval
else if (value >= values.back()) {
retval.ind_[0] = nvalues-2;
retval.ind_[1] = nvalues-1;
}
else {
//Search internal intervals
for (int i=1; i<nvalues; ++i) {
if (values[i] >= value) {
retval.ind_[0] = i-1;
retval.ind_[1] = i;
break;
}
}
}
const double start = values[retval.ind_[0]];
const double end = values[retval.ind_[1]];
//Find interpolation ratio
if (end > start) {
//FIXME: Possible source for floating point error here if value and floor are large,
//but very close to each other
retval.inv_dist_ = 1.0 / (end-start);
retval.factor_ = (value-start) * retval.inv_dist_;
}
else {
retval.inv_dist_ = 0.0;
retval.factor_ = 0.0;
}
}
return retval;
}
/**
* An "ADB-like" structure with a single value and a set of derivatives
*/
struct VFPEvaluation {
VFPEvaluation() : value(0.0), dthp(0.0), dwfr(0.0), dgfr(0.0), dalq(0.0), dflo(0.0) {};
double value;
double dthp;
double dwfr;
double dgfr;
double dalq;
double dflo;
};
inline VFPEvaluation operator+(
VFPEvaluation lhs,
const VFPEvaluation& rhs) {
lhs.value += rhs.value;
lhs.dthp += rhs.dthp;
lhs.dwfr += rhs.dwfr;
lhs.dgfr += rhs.dgfr;
lhs.dalq += rhs.dalq;
lhs.dflo += rhs.dflo;
return lhs;
}
inline VFPEvaluation operator-(
VFPEvaluation lhs,
const VFPEvaluation& rhs) {
lhs.value -= rhs.value;
lhs.dthp -= rhs.dthp;
lhs.dwfr -= rhs.dwfr;
lhs.dgfr -= rhs.dgfr;
lhs.dalq -= rhs.dalq;
lhs.dflo -= rhs.dflo;
return lhs;
}
inline VFPEvaluation operator*(
double lhs,
const VFPEvaluation& rhs) {
VFPEvaluation retval;
retval.value = rhs.value * lhs;
retval.dthp = rhs.dthp * lhs;
retval.dwfr = rhs.dwfr * lhs;
retval.dgfr = rhs.dgfr * lhs;
retval.dalq = rhs.dalq * lhs;
retval.dflo = rhs.dflo * lhs;
return retval;
}
/**
* Helper function which interpolates data using the indices etc. given in the inputs.
*/
#ifdef __GNUC__
#pragma GCC push_options
#pragma GCC optimize ("unroll-loops")
#endif
inline VFPEvaluation interpolate(
const VFPProdTable::array_type& array,
const InterpData& flo_i,
const InterpData& thp_i,
const InterpData& wfr_i,
const InterpData& gfr_i,
const InterpData& alq_i) {
//Values and derivatives in a 5D hypercube
VFPEvaluation nn[2][2][2][2][2];
//Pick out nearest neighbors (nn) to our evaluation point
//This is not really required, but performance-wise it may pay off, since the 32-elements
//we copy to (nn) will fit better in cache than the full original table for the
//interpolation below.
//The following ladder of for loops will presumably be unrolled by a reasonable compiler.
for (int t=0; t<=1; ++t) {
for (int w=0; w<=1; ++w) {
for (int g=0; g<=1; ++g) {
for (int a=0; a<=1; ++a) {
for (int f=0; f<=1; ++f) {
//Shorthands for indexing
const int ti = thp_i.ind_[t];
const int wi = wfr_i.ind_[w];
const int gi = gfr_i.ind_[g];
const int ai = alq_i.ind_[a];
const int fi = flo_i.ind_[f];
//Copy element
nn[t][w][g][a][f].value = array[ti][wi][gi][ai][fi];
}
}
}
}
}
//Calculate derivatives
//Note that the derivative of the two end points of a line aligned with the
//"axis of the derivative" are equal
for (int i=0; i<=1; ++i) {
for (int j=0; j<=1; ++j) {
for (int k=0; k<=1; ++k) {
for (int l=0; l<=1; ++l) {
nn[0][i][j][k][l].dthp = (nn[1][i][j][k][l].value - nn[0][i][j][k][l].value) * thp_i.inv_dist_;
nn[i][0][j][k][l].dwfr = (nn[i][1][j][k][l].value - nn[i][0][j][k][l].value) * wfr_i.inv_dist_;
nn[i][j][0][k][l].dgfr = (nn[i][j][1][k][l].value - nn[i][j][0][k][l].value) * gfr_i.inv_dist_;
nn[i][j][k][0][l].dalq = (nn[i][j][k][1][l].value - nn[i][j][k][0][l].value) * alq_i.inv_dist_;
nn[i][j][k][l][0].dflo = (nn[i][j][k][l][1].value - nn[i][j][k][l][0].value) * flo_i.inv_dist_;
nn[1][i][j][k][l].dthp = nn[0][i][j][k][l].dthp;
nn[i][1][j][k][l].dwfr = nn[i][0][j][k][l].dwfr;
nn[i][j][1][k][l].dgfr = nn[i][j][0][k][l].dgfr;
nn[i][j][k][1][l].dalq = nn[i][j][k][0][l].dalq;
nn[i][j][k][l][1].dflo = nn[i][j][k][l][0].dflo;
}
}
}
}
double t1, t2; //interpolation variables, so that t1 = (1-t) and t2 = t.
// Remove dimensions one by one
// Example: going from 3D to 2D to 1D, we start by interpolating along
// the z axis first, leaving a 2D problem. Then interpolating along the y
// axis, leaving a 1D, problem, etc.
t2 = flo_i.factor_;
t1 = (1.0-t2);
for (int t=0; t<=1; ++t) {
for (int w=0; w<=1; ++w) {
for (int g=0; g<=1; ++g) {
for (int a=0; a<=1; ++a) {
nn[t][w][g][a][0] = t1*nn[t][w][g][a][0] + t2*nn[t][w][g][a][1];
}
}
}
}
t2 = alq_i.factor_;
t1 = (1.0-t2);
for (int t=0; t<=1; ++t) {
for (int w=0; w<=1; ++w) {
for (int g=0; g<=1; ++g) {
nn[t][w][g][0][0] = t1*nn[t][w][g][0][0] + t2*nn[t][w][g][1][0];
}
}
}
t2 = gfr_i.factor_;
t1 = (1.0-t2);
for (int t=0; t<=1; ++t) {
for (int w=0; w<=1; ++w) {
nn[t][w][0][0][0] = t1*nn[t][w][0][0][0] + t2*nn[t][w][1][0][0];
}
}
t2 = wfr_i.factor_;
t1 = (1.0-t2);
for (int t=0; t<=1; ++t) {
nn[t][0][0][0][0] = t1*nn[t][0][0][0][0] + t2*nn[t][1][0][0][0];
}
t2 = thp_i.factor_;
t1 = (1.0-t2);
nn[0][0][0][0][0] = t1*nn[0][0][0][0][0] + t2*nn[1][0][0][0][0];
return nn[0][0][0][0][0];
}
/**
* This basically models interpolate(VFPProdTable::array_type, ...)
* which performs 5D interpolation, but here for the 2D case only
*/
inline VFPEvaluation interpolate(
const VFPInjTable::array_type& array,
const InterpData& flo_i,
const InterpData& thp_i) {
//Values and derivatives in a 2D plane
VFPEvaluation nn[2][2];
//Pick out nearest neighbors (nn) to our evaluation point
//The following ladder of for loops will presumably be unrolled by a reasonable compiler.
for (int t=0; t<=1; ++t) {
for (int f=0; f<=1; ++f) {
//Shorthands for indexing
const int ti = thp_i.ind_[t];
const int fi = flo_i.ind_[f];
//Copy element
nn[t][f].value = array[ti][fi];
}
}
//Calculate derivatives
//Note that the derivative of the two end points of a line aligned with the
//"axis of the derivative" are equal
for (int i=0; i<=1; ++i) {
nn[0][i].dthp = (nn[1][i].value - nn[0][i].value) * thp_i.inv_dist_;
nn[i][0].dwfr = -1e100;
nn[i][0].dgfr = -1e100;
nn[i][0].dalq = -1e100;
nn[i][0].dflo = (nn[i][1].value - nn[i][0].value) * flo_i.inv_dist_;
nn[1][i].dthp = nn[0][i].dthp;
nn[i][1].dwfr = nn[i][0].dwfr;
nn[i][1].dgfr = nn[i][0].dgfr;
nn[i][1].dalq = nn[i][0].dalq;
nn[i][1].dflo = nn[i][0].dflo;
}
double t1, t2; //interpolation variables, so that t1 = (1-t) and t2 = t.
// Go from 2D to 1D
t2 = flo_i.factor_;
t1 = (1.0-t2);
nn[0][0] = t1*nn[0][0] + t2*nn[0][1];
nn[1][0] = t1*nn[1][0] + t2*nn[1][1];
// Go from line to point on line
t2 = thp_i.factor_;
t1 = (1.0-t2);
nn[0][0] = t1*nn[0][0] + t2*nn[1][0];
return nn[0][0];
}
#ifdef __GNUC__
#pragma GCC pop_options //unroll loops
#endif
inline VFPEvaluation bhp(const VFPProdTable* table,
const double& aqua,
const double& liquid,
const double& vapour,
const double& thp,
const double& alq) {
//Find interpolation variables
double flo = detail::getFlo(aqua, liquid, vapour, table->getFloType());
double wfr = detail::getWFR(aqua, liquid, vapour, table->getWFRType());
double gfr = detail::getGFR(aqua, liquid, vapour, table->getGFRType());
//First, find the values to interpolate between
//Recall that flo is negative in Opm, so switch sign.
auto flo_i = detail::findInterpData(-flo, table->getFloAxis());
auto thp_i = detail::findInterpData( thp, table->getTHPAxis());
auto wfr_i = detail::findInterpData( wfr, table->getWFRAxis());
auto gfr_i = detail::findInterpData( gfr, table->getGFRAxis());
auto alq_i = detail::findInterpData( alq, table->getALQAxis());
detail::VFPEvaluation retval = detail::interpolate(table->getTable(), flo_i, thp_i, wfr_i, gfr_i, alq_i);
return retval;
}
inline VFPEvaluation bhp(const VFPInjTable* table,
const double& aqua,
const double& liquid,
const double& vapour,
const double& thp) {
//Find interpolation variables
double flo = detail::getFlo(aqua, liquid, vapour, table->getFloType());
//First, find the values to interpolate between
auto flo_i = detail::findInterpData(flo, table->getFloAxis());
auto thp_i = detail::findInterpData(thp, table->getTHPAxis());
//Then perform the interpolation itself
detail::VFPEvaluation retval = detail::interpolate(table->getTable(), flo_i, thp_i);
return retval;
}
/**
* Returns the table from the map if found, or throws an exception
*/
template <typename T>
const T* getTable(const std::map<int, T*> tables, int table_id) {
auto entry = tables.find(table_id);
if (entry == tables.end()) {
OPM_THROW(std::invalid_argument, "Nonexistent table " << table_id << " referenced.");
}
else {
return entry->second;
}
}
/**
* Sets block_pattern to be the "union of x.blockPattern() and block_pattern".
*/
inline void extendBlockPattern(const ADB& x, std::vector<int>& block_pattern) {
std::vector<int> x_block_pattern = x.blockPattern();
if (x_block_pattern.empty()) {
return;
}
else {
if (block_pattern.empty()) {
block_pattern = x_block_pattern;
return;
}
else {
if (x_block_pattern != block_pattern) {
OPM_THROW(std::logic_error, "Block patterns do not match");
}
}
}
}
/**
* Finds the common block pattern for all inputs
*/
inline std::vector<int> commonBlockPattern(
const ADB& x1,
const ADB& x2,
const ADB& x3,
const ADB& x4) {
std::vector<int> block_pattern;
extendBlockPattern(x1, block_pattern);
extendBlockPattern(x2, block_pattern);
extendBlockPattern(x3, block_pattern);
extendBlockPattern(x4, block_pattern);
return block_pattern;
}
inline std::vector<int> commonBlockPattern(
const ADB& x1,
const ADB& x2,
const ADB& x3,
const ADB& x4,
const ADB& x5) {
std::vector<int> block_pattern = commonBlockPattern(x1, x2, x3, x4);
extendBlockPattern(x5, block_pattern);
return block_pattern;
}
/**
* Returns the type variable for FLO/GFR/WFR for production tables
*/
template <typename TYPE, typename TABLE>
TYPE getType(const TABLE* table);
template <>
inline
VFPProdTable::FLO_TYPE getType(const VFPProdTable* table) {
return table->getFloType();
}
template <>
inline
VFPProdTable::WFR_TYPE getType(const VFPProdTable* table) {
return table->getWFRType();
}
template <>
inline
VFPProdTable::GFR_TYPE getType(const VFPProdTable* table) {
return table->getGFRType();
}
/**
* Returns the type variable for FLO for injection tables
*/
template <>
inline
VFPInjTable::FLO_TYPE getType(const VFPInjTable* table) {
return table->getFloType();
}
/**
* Returns the actual ADB for the type of FLO/GFR/WFR type
*/
template <typename TYPE>
ADB getValue(
const ADB& aqua,
const ADB& liquid,
const ADB& vapour, TYPE type);
template <>
inline
ADB getValue(
const ADB& aqua,
const ADB& liquid,
const ADB& vapour,
VFPProdTable::FLO_TYPE type) {
return detail::getFlo(aqua, liquid, vapour, type);
}
template <>
inline
ADB getValue(
const ADB& aqua,
const ADB& liquid,
const ADB& vapour,
VFPProdTable::WFR_TYPE type) {
return detail::getWFR(aqua, liquid, vapour, type);
}
template <>
inline
ADB getValue(
const ADB& aqua,
const ADB& liquid,
const ADB& vapour,
VFPProdTable::GFR_TYPE type) {
return detail::getGFR(aqua, liquid, vapour, type);
}
template <>
inline
ADB getValue(
const ADB& aqua,
const ADB& liquid,
const ADB& vapour,
VFPInjTable::FLO_TYPE type) {
return detail::getFlo(aqua, liquid, vapour, type);
}
/**
* Given m wells and n types of VFP variables (e.g., FLO = {FLO_OIL, FLO_LIQ}
* this function combines the n types of ADB objects, so that each of the
* m wells gets the right ADB.
* @param TYPE Type of variable to return, e.g., FLO_TYPE, WFR_TYPE, GFR_TYPE
* @param TABLE Type of table to use, e.g., VFPInjTable, VFPProdTable.
*/
template <typename TYPE, typename TABLE>
ADB combineADBVars(const std::vector<const TABLE*>& well_tables,
const ADB& aqua,
const ADB& liquid,
const ADB& vapour) {
const int num_wells = static_cast<int>(well_tables.size());
assert(aqua.size() == num_wells);
assert(liquid.size() == num_wells);
assert(vapour.size() == num_wells);
//Caching variable for flo/wfr/gfr
std::map<TYPE, ADB> map;
//Indexing variable used when combining the different ADB types
std::map<TYPE, std::vector<int> > elems;
//Compute all of the different ADB types,
//and record which wells use which types
for (int i=0; i<num_wells; ++i) {
const TABLE* table = well_tables[i];
//Only do something if this well is under THP control
if (table != NULL) {
TYPE type = getType<TYPE>(table);
//"Caching" of flo_type etc: Only calculate used types
//Create type if it does not exist
if (map.find(type) == map.end()) {
map.insert(std::pair<TYPE, ADB>(
type,
detail::getValue<TYPE>(aqua, liquid, vapour, type)
));
}
//Add the index for assembly later in gather_vars
elems[type].push_back(i);
}
}
//Loop over all types of ADB variables, and combine them
//so that each well gets the proper variable
ADB retval = ADB::constant(ADB::V::Zero(num_wells));
for (const auto& entry : elems) {
const auto& key = entry.first;
const auto& value = entry.second;
//Get the ADB for this type of variable
assert(map.find(key) != map.end());
const ADB& values = map.find(key)->second;
//Get indices to all elements that should use this ADB
const std::vector<int>& current = value;
//Add these elements to retval
retval = retval + superset(subset(values, current), current, values.size());
}
return retval;
}
/**
* Helper function that finds x for a given value of y for a line
* *NOTE ORDER OF ARGUMENTS*
*/
inline double findX(const double& x0,
const double& x1,
const double& y0,
const double& y1,
const double& y) {
const double dx = x1 - x0;
const double dy = y1 - y0;
/**
* y = y0 + (dy / dx) * (x - x0)
* => x = x0 + (y - y0) * (dx / dy)
*
* If dy is zero, use x1 as the value.
*/
double x = 0.0;
if (dy != 0.0) {
x = x0 + (y-y0) * (dx/dy);
}
else {
x = x1;
}
return x;
}
/**
* This function finds the value of THP given a specific BHP.
* Essentially:
* Given the function f(thp_array(x)) = bhp_array(x), which is piecewise linear,
* find thp so that f(thp) = bhp.
*/
inline double findTHP(
const std::vector<double>& bhp_array,
const std::vector<double>& thp_array,
double bhp) {
int nthp = thp_array.size();
double thp = -1e100;
//Check that our thp axis is sorted
assert(std::is_sorted(thp_array.begin(), thp_array.end()));
/**
* Our *interpolated* bhp_array will be montonic increasing for increasing
* THP if our input BHP values are monotonic increasing for increasing
* THP values. However, if we have to *extrapolate* along any of the other
* axes, this guarantee holds no more, and bhp_array may be "random"
*/
if (std::is_sorted(bhp_array.begin(), bhp_array.end())) {
//Target bhp less than all values in array, extrapolate
if (bhp <= bhp_array[0]) {
//TODO: LOG extrapolation
const double& x0 = thp_array[0];
const double& x1 = thp_array[1];
const double& y0 = bhp_array[0];
const double& y1 = bhp_array[1];
thp = detail::findX(x0, x1, y0, y1, bhp);
}
//Target bhp greater than all values in array, extrapolate
else if (bhp > bhp_array[nthp-1]) {
//TODO: LOG extrapolation
const double& x0 = thp_array[nthp-2];
const double& x1 = thp_array[nthp-1];
const double& y0 = bhp_array[nthp-2];
const double& y1 = bhp_array[nthp-1];
thp = detail::findX(x0, x1, y0, y1, bhp);
}
//Target bhp within table ranges, interpolate
else {
//Loop over the values and find min(bhp_array(thp)) == bhp
//so that we maximize the rate.
//Find i so that bhp_array[i-1] <= bhp <= bhp_array[i];
//Assuming a small number of values in bhp_array, this should be quite
//efficient. Other strategies might be bisection, etc.
int i=0;
bool found = false;
for (; i<nthp-1; ++i) {
const double& y0 = bhp_array[i ];
const double& y1 = bhp_array[i+1];
if (y0 < bhp && bhp <= y1) {
found = true;
break;
}
}
//Canary in a coal mine: shouldn't really be required
assert(found == true);
const double& x0 = thp_array[i ];
const double& x1 = thp_array[i+1];
const double& y0 = bhp_array[i ];
const double& y1 = bhp_array[i+1];
thp = detail::findX(x0, x1, y0, y1, bhp);
}
}
//bhp_array not sorted, raw search.
else {
//Find i so that bhp_array[i-1] <= bhp <= bhp_array[i];
//Since the BHP values might not be sorted, first search within
//our interpolation values, and then try to extrapolate.
int i=0;
bool found = false;
for (; i<nthp-1; ++i) {
const double& y0 = bhp_array[i ];
const double& y1 = bhp_array[i+1];
if (y0 < bhp && bhp <= y1) {
found = true;
break;
}
}
if (found) {
const double& x0 = thp_array[i ];
const double& x1 = thp_array[i+1];
const double& y0 = bhp_array[i ];
const double& y1 = bhp_array[i+1];
thp = detail::findX(x0, x1, y0, y1, bhp);
}
else if (bhp <= bhp_array[0]) {
//TODO: LOG extrapolation
const double& x0 = thp_array[0];
const double& x1 = thp_array[1];
const double& y0 = bhp_array[0];
const double& y1 = bhp_array[1];
thp = detail::findX(x0, x1, y0, y1, bhp);
}
//Target bhp greater than all values in array, extrapolate
else if (bhp > bhp_array[nthp-1]) {
//TODO: LOG extrapolation
const double& x0 = thp_array[nthp-2];
const double& x1 = thp_array[nthp-1];
const double& y0 = bhp_array[nthp-2];
const double& y1 = bhp_array[nthp-1];
thp = detail::findX(x0, x1, y0, y1, bhp);
}
else {
OPM_THROW(std::logic_error, "Programmer error: Unable to find THP in THP array");
}
}
return thp;
}
} // namespace detail
} // namespace
#endif /* OPM_AUTODIFF_VFPHELPERS_HPP_ */

238
opm/autodiff/VFPInjProperties.cpp Normal file
View File

@@ -0,0 +1,238 @@
/*
Copyright 2015 SINTEF ICT, Applied Mathematics.
This file is part of the Open Porous Media project (OPM).
OPM is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OPM is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#include "config.h"
#include <opm/autodiff/VFPInjProperties.hpp>
#include <opm/parser/eclipse/EclipseState/Tables/VFPProdTable.hpp>
#include <opm/core/props/BlackoilPhases.hpp>
#include <opm/core/utility/ErrorMacros.hpp>
#include <opm/autodiff/AutoDiffHelpers.hpp>
#include <opm/autodiff/VFPHelpers.hpp>
namespace Opm {
namespace detail {
} //Namespace detail
VFPInjProperties::VFPInjProperties() {
}
VFPInjProperties::VFPInjProperties(const VFPInjTable* table){
m_tables[table->getTableNum()] = table;
}
VFPInjProperties::VFPInjProperties(const std::map<int, VFPInjTable>& tables) {
for (const auto& table : tables) {
m_tables[table.first] = &table.second;
}
}
VFPInjProperties::ADB VFPInjProperties::bhp(const std::vector<int>& table_id,
const Wells& wells,
const ADB& qs,
const ADB& thp) const {
const int np = wells.number_of_phases;
const int nw = wells.number_of_wells;
//Short-hands for water / oil / gas phases
//TODO enable support for two-phase.
assert(np == 3);
const ADB& w = subset(qs, Span(nw, 1, BlackoilPhases::Aqua*nw));
const ADB& o = subset(qs, Span(nw, 1, BlackoilPhases::Liquid*nw));
const ADB& g = subset(qs, Span(nw, 1, BlackoilPhases::Vapour*nw));
return bhp(table_id, w, o, g, thp);
}
VFPInjProperties::ADB VFPInjProperties::bhp(const std::vector<int>& table_id,
const ADB& aqua,
const ADB& liquid,
const ADB& vapour,
const ADB& thp) const {
const int nw = thp.size();
std::vector<int> block_pattern = detail::commonBlockPattern(aqua, liquid, vapour, thp);
assert(static_cast<int>(table_id.size()) == nw);
assert(aqua.size() == nw);
assert(liquid.size() == nw);
assert(vapour.size() == nw);
assert(thp.size() == nw);
//Allocate data for bhp's and partial derivatives
ADB::V value = ADB::V::Zero(nw);
ADB::V dthp = ADB::V::Zero(nw);
ADB::V dflo = ADB::V::Zero(nw);
//Get the table for each well
std::vector<const VFPInjTable*> well_tables(nw, nullptr);
for (int i=0; i<nw; ++i) {
if (table_id[i] > 0) {
well_tables[i] = detail::getTable(m_tables, table_id[i]);
}
}
//Get the right FLO variable for each well as a single ADB
const ADB flo = detail::combineADBVars<VFPInjTable::FLO_TYPE>(well_tables, aqua, liquid, vapour);
//Compute the BHP for each well independently
for (int i=0; i<nw; ++i) {
const VFPInjTable* table = well_tables[i];
if (table != nullptr) {
//First, find the values to interpolate between
auto flo_i = detail::findInterpData(flo.value()[i], table->getFloAxis());
auto thp_i = detail::findInterpData(thp.value()[i], table->getTHPAxis());
detail::VFPEvaluation bhp_val = detail::interpolate(table->getTable(), flo_i, thp_i);
value[i] = bhp_val.value;
dthp[i] = bhp_val.dthp;
dflo[i] = bhp_val.dflo;
}
else {
value[i] = -1e100; //Signal that this value has not been calculated properly, due to "missing" table
}
}
//Create diagonal matrices from ADB::Vs
ADB::M dthp_diag = spdiag(dthp);
ADB::M dflo_diag = spdiag(dflo);
//Calculate the Jacobians
const int num_blocks = block_pattern.size();
std::vector<ADB::M> jacs(num_blocks);
for (int block = 0; block < num_blocks; ++block) {
//Could have used fastSparseProduct and temporary variables
//but may not save too much on that.
jacs[block] = ADB::M(nw, block_pattern[block]);
if (!thp.derivative().empty()) {
jacs[block] += dthp_diag * thp.derivative()[block];
}
if (!flo.derivative().empty()) {
jacs[block] += dflo_diag * flo.derivative()[block];
}
}
ADB retval = ADB::function(std::move(value), std::move(jacs));
return retval;
}
double VFPInjProperties::bhp(int table_id,
const double& aqua,
const double& liquid,
const double& vapour,
const double& thp) const {
const VFPInjTable* table = detail::getTable(m_tables, table_id);
detail::VFPEvaluation retval = detail::bhp(table, aqua, liquid, vapour, thp);
return retval.value;
}
double VFPInjProperties::thp(int table_id,
const double& aqua,
const double& liquid,
const double& vapour,
const double& bhp) const {
const VFPInjTable* table = detail::getTable(m_tables, table_id);
const VFPInjTable::array_type& data = table->getTable();
//Find interpolation variables
double flo = detail::getFlo(aqua, liquid, vapour, table->getFloType());
const std::vector<double> thp_array = table->getTHPAxis();
int nthp = thp_array.size();
/**
* Find the function bhp_array(thp) by creating a 1D view of the data
* by interpolating for every value of thp. This might be somewhat
* expensive, but let us assome that nthp is small
*/
auto flo_i = detail::findInterpData(flo, table->getFloAxis());
std::vector<double> bhp_array(nthp);
for (int i=0; i<nthp; ++i) {
auto thp_i = detail::findInterpData(thp_array[i], thp_array);
bhp_array[i] = detail::interpolate(data, flo_i, thp_i).value;
}
double thp = detail::findTHP(bhp_array, thp_array, bhp);
return thp;
}
const VFPInjTable* VFPInjProperties::getTable(const int table_id) const {
return detail::getTable(m_tables, table_id);
}
} //Namespace Opm

153
opm/autodiff/VFPInjProperties.hpp Normal file
View File

@@ -0,0 +1,153 @@
/*
Copyright 2015 SINTEF ICT, Applied Mathematics.
This file is part of the Open Porous Media project (OPM).
OPM is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OPM is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef OPM_AUTODIFF_VFPINJPROPERTIES_HPP_
#define OPM_AUTODIFF_VFPINJPROPERTIES_HPP_
#include <opm/parser/eclipse/EclipseState/Tables/VFPInjTable.hpp>
#include <opm/core/wells.h>
#include <opm/autodiff/AutoDiffBlock.hpp>
#include <vector>
#include <map>
namespace Opm {
class VFPInjProperties {
public:
typedef AutoDiffBlock<double> ADB;
/**
* Empty constructor
*/
VFPInjProperties();
/**
* Constructor
* Takes *no* ownership of data.
* @param inj_table A *single* VFPINJ table
*/
explicit VFPInjProperties(const VFPInjTable* inj_table);
/**
* Constructor
* Takes *no* ownership of data.
* @param inj_tables A map of different VFPINJ tables.
*/
explicit VFPInjProperties(const std::map<int, VFPInjTable>& inj_tables);
/**
* Linear interpolation of bhp as function of the input parameters.
* @param table_id Table number to use
* @param wells Wells structure with information about wells in qs
* @param qs Flow quantities
* @param thp Tubing head pressure
*
* @return The bottom hole pressure, interpolated/extrapolated linearly using
* the above parameters from the values in the input table.
*/
ADB bhp(const std::vector<int>& table_id,
const Wells& wells,
const ADB& qs,
const ADB& thp) const;
/**
* Linear interpolation of bhp as a function of the input parameters given as ADBs
* Each entry corresponds typically to one well.
* @param table_id Table number to use. A negative entry (e.g., -1)
* will indicate that no table is used, and the corresponding
* BHP will be calculated as a constant -1e100.
* @param aqua Water phase
* @param liquid Oil phase
* @param vapour Gas phase
* @param thp Tubing head pressure
*
* @return The bottom hole pressure, interpolated/extrapolated linearly using
* the above parameters from the values in the input table, for each entry in the
* input ADB objects.
*/
ADB bhp(const std::vector<int>& table_id,
const ADB& aqua,
const ADB& liquid,
const ADB& vapour,
const ADB& thp) const;
/**
* Linear interpolation of bhp as a function of the input parameters
* @param table_id Table number to use
* @param aqua Water phase
* @param liquid Oil phase
* @param vapour Gas phase
* @param thp Tubing head pressure
*
* @return The bottom hole pressure, interpolated/extrapolated linearly using
* the above parameters from the values in the input table.
*/
double bhp(int table_id,
const double& aqua,
const double& liquid,
const double& vapour,
const double& thp) const;
/**
* Linear interpolation of thp as a function of the input parameters
* @param table_id Table number to use
* @param aqua Water phase
* @param liquid Oil phase
* @param vapour Gas phase
* @param bhp Bottom hole pressure
*
* @return The tubing hole pressure, interpolated/extrapolated linearly using
* the above parameters from the values in the input table.
*/
double thp(int table_id,
const double& aqua,
const double& liquid,
const double& vapour,
const double& bhp) const;
/**
* Returns the table associated with the ID, or throws an exception if
* the table does not exist
*/
const VFPInjTable* getTable(const int table_id) const;
/**
* Returns true if no vfp tables are in the current map
*/
bool empty() const {
return m_tables.empty();
}
private:
// Map which connects the table number with the table itself
std::map<int, const VFPInjTable*> m_tables;
};
} //namespace
#endif /* OPM_AUTODIFF_VFPINJPROPERTIES_HPP_ */

250
opm/autodiff/VFPProdProperties.cpp Normal file
View File

@@ -0,0 +1,250 @@
/*
Copyright 2015 SINTEF ICT, Applied Mathematics.
This file is part of the Open Porous Media project (OPM).
OPM is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OPM is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#include "config.h"
#include <opm/autodiff/VFPProdProperties.hpp>
#include <opm/parser/eclipse/EclipseState/Tables/VFPProdTable.hpp>
#include <opm/core/props/BlackoilPhases.hpp>
#include <opm/core/utility/ErrorMacros.hpp>
#include <opm/autodiff/AutoDiffHelpers.hpp>
#include <opm/autodiff/VFPHelpers.hpp>
namespace Opm {
VFPProdProperties::VFPProdProperties() {
}
VFPProdProperties::VFPProdProperties(const VFPProdTable* table){
m_tables[table->getTableNum()] = table;
}
VFPProdProperties::VFPProdProperties(const std::map<int, VFPProdTable>& tables) {
for (const auto& table : tables) {
m_tables[table.first] = &table.second;
}
}
VFPProdProperties::ADB VFPProdProperties::bhp(const std::vector<int>& table_id,
const Wells& wells,
const ADB& qs,
const ADB& thp,
const ADB& alq) const {
const int np = wells.number_of_phases;
const int nw = wells.number_of_wells;
//Short-hands for water / oil / gas phases
//TODO enable support for two-phase.
assert(np == 3);
const ADB& w = subset(qs, Span(nw, 1, BlackoilPhases::Aqua*nw));
const ADB& o = subset(qs, Span(nw, 1, BlackoilPhases::Liquid*nw));
const ADB& g = subset(qs, Span(nw, 1, BlackoilPhases::Vapour*nw));
return bhp(table_id, w, o, g, thp, alq);
}
VFPProdProperties::ADB VFPProdProperties::bhp(const std::vector<int>& table_id,
const ADB& aqua,
const ADB& liquid,
const ADB& vapour,
const ADB& thp,
const ADB& alq) const {
const int nw = thp.size();
std::vector<int> block_pattern = detail::commonBlockPattern(aqua, liquid, vapour, thp, alq);
assert(static_cast<int>(table_id.size()) == nw);
assert(aqua.size() == nw);
assert(liquid.size() == nw);
assert(vapour.size() == nw);
assert(thp.size() == nw);
assert(alq.size() == nw);
//Allocate data for bhp's and partial derivatives
ADB::V value = ADB::V::Zero(nw);
ADB::V dthp = ADB::V::Zero(nw);
ADB::V dwfr = ADB::V::Zero(nw);
ADB::V dgfr = ADB::V::Zero(nw);
ADB::V dalq = ADB::V::Zero(nw);
ADB::V dflo = ADB::V::Zero(nw);
//Get the table for each well
std::vector<const VFPProdTable*> well_tables(nw, nullptr);
for (int i=0; i<nw; ++i) {
if (table_id[i] >= 0) {
well_tables[i] = detail::getTable(m_tables, table_id[i]);
}
}
//Get the right FLO/GFR/WFR variable for each well as a single ADB
const ADB flo = detail::combineADBVars<VFPProdTable::FLO_TYPE>(well_tables, aqua, liquid, vapour);
const ADB wfr = detail::combineADBVars<VFPProdTable::WFR_TYPE>(well_tables, aqua, liquid, vapour);
const ADB gfr = detail::combineADBVars<VFPProdTable::GFR_TYPE>(well_tables, aqua, liquid, vapour);
//Compute the BHP for each well independently
for (int i=0; i<nw; ++i) {
const VFPProdTable* table = well_tables[i];
if (table != nullptr) {
//First, find the values to interpolate between
//Value of FLO is negative in OPM for producers, but positive in VFP table
auto flo_i = detail::findInterpData(-flo.value()[i], table->getFloAxis());
auto thp_i = detail::findInterpData( thp.value()[i], table->getTHPAxis());
auto wfr_i = detail::findInterpData( wfr.value()[i], table->getWFRAxis());
auto gfr_i = detail::findInterpData( gfr.value()[i], table->getGFRAxis());
auto alq_i = detail::findInterpData( alq.value()[i], table->getALQAxis());
detail::VFPEvaluation bhp_val = detail::interpolate(table->getTable(), flo_i, thp_i, wfr_i, gfr_i, alq_i);
value[i] = bhp_val.value;
dthp[i] = bhp_val.dthp;
dwfr[i] = bhp_val.dwfr;
dgfr[i] = bhp_val.dgfr;
dalq[i] = bhp_val.dalq;
dflo[i] = bhp_val.dflo;
}
else {
value[i] = -1e100; //Signal that this value has not been calculated properly, due to "missing" table
}
}
//Create diagonal matrices from ADB::Vs
ADB::M dthp_diag = spdiag(dthp);
ADB::M dwfr_diag = spdiag(dwfr);
ADB::M dgfr_diag = spdiag(dgfr);
ADB::M dalq_diag = spdiag(dalq);
ADB::M dflo_diag = spdiag(dflo);
//Calculate the Jacobians
const int num_blocks = block_pattern.size();
std::vector<ADB::M> jacs(num_blocks);
for (int block = 0; block < num_blocks; ++block) {
//Could have used fastSparseProduct and temporary variables
//but may not save too much on that.
jacs[block] = ADB::M(nw, block_pattern[block]);
if (!thp.derivative().empty()) {
jacs[block] += dthp_diag * thp.derivative()[block];
}
if (!wfr.derivative().empty()) {
jacs[block] += dwfr_diag * wfr.derivative()[block];
}
if (!gfr.derivative().empty()) {
jacs[block] += dgfr_diag * gfr.derivative()[block];
}
if (!alq.derivative().empty()) {
jacs[block] += dalq_diag * alq.derivative()[block];
}
if (!flo.derivative().empty()) {
jacs[block] -= dflo_diag * flo.derivative()[block];
}
}
ADB retval = ADB::function(std::move(value), std::move(jacs));
return retval;
}
double VFPProdProperties::bhp(int table_id,
const double& aqua,
const double& liquid,
const double& vapour,
const double& thp,
const double& alq) const {
const VFPProdTable* table = detail::getTable(m_tables, table_id);
detail::VFPEvaluation retval = detail::bhp(table, aqua, liquid, vapour, thp, alq);
return retval.value;
}
double VFPProdProperties::thp(int table_id,
const double& aqua,
const double& liquid,
const double& vapour,
const double& bhp,
const double& alq) const {
const VFPProdTable* table = detail::getTable(m_tables, table_id);
const VFPProdTable::array_type& data = table->getTable();
//Find interpolation variables
double flo = detail::getFlo(aqua, liquid, vapour, table->getFloType());
double wfr = detail::getWFR(aqua, liquid, vapour, table->getWFRType());
double gfr = detail::getGFR(aqua, liquid, vapour, table->getGFRType());
const std::vector<double> thp_array = table->getTHPAxis();
int nthp = thp_array.size();
/**
* Find the function bhp_array(thp) by creating a 1D view of the data
* by interpolating for every value of thp. This might be somewhat
* expensive, but let us assome that nthp is small
* Recall that flo is negative in Opm, so switch the sign
*/
auto flo_i = detail::findInterpData(-flo, table->getFloAxis());
auto wfr_i = detail::findInterpData( wfr, table->getWFRAxis());
auto gfr_i = detail::findInterpData( gfr, table->getGFRAxis());
auto alq_i = detail::findInterpData( alq, table->getALQAxis());
std::vector<double> bhp_array(nthp);
for (int i=0; i<nthp; ++i) {
auto thp_i = detail::findInterpData(thp_array[i], thp_array);
bhp_array[i] = detail::interpolate(data, flo_i, thp_i, wfr_i, gfr_i, alq_i).value;
}
double thp = detail::findTHP(bhp_array, thp_array, bhp);
return thp;
}
const VFPProdTable* VFPProdProperties::getTable(const int table_id) const {
return detail::getTable(m_tables, table_id);
}
}

164
opm/autodiff/VFPProdProperties.hpp Normal file
View File

@@ -0,0 +1,164 @@
/*
Copyright 2015 SINTEF ICT, Applied Mathematics.
This file is part of the Open Porous Media project (OPM).
OPM is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OPM is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef OPM_AUTODIFF_VFPPRODPROPERTIES_HPP_
#define OPM_AUTODIFF_VFPPRODPROPERTIES_HPP_
#include <opm/parser/eclipse/EclipseState/Tables/VFPProdTable.hpp>
#include <opm/core/wells.h>
#include <opm/autodiff/AutoDiffBlock.hpp>
#include <vector>
#include <map>
namespace Opm {
/**
* Class which linearly interpolates BHP as a function of rate, tubing head pressure,
* water fraction, gas fraction, and artificial lift for production VFP tables, and similarly
* the BHP as a function of the rate and tubing head pressure.
*/
class VFPProdProperties {
public:
typedef AutoDiffBlock<double> ADB;
/**
* Empty constructor
*/
VFPProdProperties();
/**
* Constructor
* Takes *no* ownership of data.
* @param prod_table A *single* VFPPROD table
*/
explicit VFPProdProperties(const VFPProdTable* prod_table);
/**
* Constructor
* Takes *no* ownership of data.
* @param prod_tables A map of different VFPPROD tables.
*/
explicit VFPProdProperties(const std::map<int, VFPProdTable>& prod_tables);
/**
* Linear interpolation of bhp as function of the input parameters.
* @param table_id Table number to use
* @param wells Wells structure with information about wells in qs
* @param qs Flow quantities
* @param thp Tubing head pressure
* @param alq Artificial lift or other parameter
*
* @return The bottom hole pressure, interpolated/extrapolated linearly using
* the above parameters from the values in the input table.
*/
ADB bhp(const std::vector<int>& table_id,
const Wells& wells,
const ADB& qs,
const ADB& thp,
const ADB& alq) const;
/**
* Linear interpolation of bhp as a function of the input parameters given as ADBs
* Each entry corresponds typically to one well.
* @param table_id Table number to use. A negative entry (e.g., -1)
* will indicate that no table is used, and the corresponding
* BHP will be calculated as a constant -1e100.
* @param aqua Water phase
* @param liquid Oil phase
* @param vapour Gas phase
* @param thp Tubing head pressure
* @param alq Artificial lift or other parameter
*
* @return The bottom hole pressure, interpolated/extrapolated linearly using
* the above parameters from the values in the input table, for each entry in the
* input ADB objects.
*/
ADB bhp(const std::vector<int>& table_id,
const ADB& aqua,
const ADB& liquid,
const ADB& vapour,
const ADB& thp,
const ADB& alq) const;
/**
* Linear interpolation of bhp as a function of the input parameters
* @param table_id Table number to use
* @param aqua Water phase
* @param liquid Oil phase
* @param vapour Gas phase
* @param thp Tubing head pressure
* @param alq Artificial lift or other parameter
*
* @return The bottom hole pressure, interpolated/extrapolated linearly using
* the above parameters from the values in the input table.
*/
double bhp(int table_id,
const double& aqua,
const double& liquid,
const double& vapour,
const double& thp,
const double& alq) const;
/**
* Linear interpolation of thp as a function of the input parameters
* @param table_id Table number to use
* @param aqua Water phase
* @param liquid Oil phase
* @param vapour Gas phase
* @param bhp Bottom hole pressure
* @param alq Artificial lift or other parameter
*
* @return The tubing hole pressure, interpolated/extrapolated linearly using
* the above parameters from the values in the input table.
*/
double thp(int table_id,
const double& aqua,
const double& liquid,
const double& vapour,
const double& bhp,
const double& alq) const;
/**
* Returns the table associated with the ID, or throws an exception if
* the table does not exist
*/
const VFPProdTable* getTable(const int table_id) const;
/**
* Returns true if no vfp tables are in the current map
*/
bool empty() const {
return m_tables.empty();
}
private:
// Map which connects the table number with the table itself
std::map<int, const VFPProdTable*> m_tables;
};
} //namespace
#endif /* OPM_AUTODIFF_VFPPRODPROPERTIES_HPP_ */

49
opm/autodiff/VFPProperties.cpp Normal file
View File

@@ -0,0 +1,49 @@
/*
Copyright 2015 SINTEF ICT, Applied Mathematics.
This file is part of the Open Porous Media project (OPM).
OPM is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OPM is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#include "config.h"
#include <opm/autodiff/VFPProperties.hpp>
#include <opm/autodiff/VFPProdProperties.hpp>
#include <opm/autodiff/VFPInjProperties.hpp>
namespace Opm {
VFPProperties::VFPProperties() {
}
VFPProperties::VFPProperties(const VFPInjTable* inj_table, const VFPProdTable* prod_table) {
if (inj_table != NULL) {
m_inj.reset(new VFPInjProperties(inj_table));
}
if (prod_table != NULL) {
m_prod.reset(new VFPProdProperties(prod_table));
}
}
VFPProperties::VFPProperties(const std::map<int, VFPInjTable>& inj_tables,
const std::map<int, VFPProdTable>& prod_tables) {
m_inj.reset(new VFPInjProperties(inj_tables));
m_prod.reset(new VFPProdProperties(prod_tables));
}
} //Namespace

80
opm/autodiff/VFPProperties.hpp Normal file
View File

@@ -0,0 +1,80 @@
/*
Copyright 2015 SINTEF ICT, Applied Mathematics.
This file is part of the Open Porous Media project (OPM).
OPM is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OPM is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef OPM_AUTODIFF_VFPPROPERTIES_HPP_
#define OPM_AUTODIFF_VFPPROPERTIES_HPP_
#include <opm/parser/eclipse/EclipseState/Tables/VFPInjTable.hpp>
#include <opm/parser/eclipse/EclipseState/Tables/VFPProdTable.hpp>
#include <map>
namespace Opm {
class VFPProdProperties;
class VFPInjProperties;
/**
* A thin wrapper class that holds one VFPProdProperties and one
* VFPInjProperties object.
*/
class VFPProperties {
public:
VFPProperties();
/**
* Constructor
* Takes *no* ownership of data.
* @param inj_table A *single* VFPINJ table or NULL (no table)
* @param prod_table A *single* VFPPROD table or NULL (no table)
*/
explicit VFPProperties(const VFPInjTable* inj_table, const VFPProdTable* prod_table);
/**
* Constructor
* Takes *no* ownership of data.
* @param inj_tables A map of different VFPINJ tables.
* @param prod_tables A map of different VFPPROD tables.
*/
VFPProperties(const std::map<int, VFPInjTable>& inj_tables,
const std::map<int, VFPProdTable>& prod_tables);
/**
* Returns the VFP properties for injection wells
*/
const VFPInjProperties* getInj() const {
return m_inj.get();
}
/**
* Returns the VFP properties for production wells
*/
const VFPProdProperties* getProd() const {
return m_prod.get();
}
private:
std::shared_ptr<VFPInjProperties> m_inj;
std::shared_ptr<VFPProdProperties> m_prod;
};
} //Namespace
#endif /* OPM_AUTODIFF_VFPPROPERTIES_HPP_ */

66
opm/autodiff/WellDensitySegmented.cpp
View File

@@ -26,16 +26,17 @@
#include <cmath>
std::vector<double>
Opm::WellDensitySegmented::computeConnectionPressureDelta(const Wells& wells,
const WellStateFullyImplicitBlackoil& wstate,
const PhaseUsage& phase_usage,
const std::vector<double>& b_perf,
const std::vector<double>& rsmax_perf,
const std::vector<double>& rvmax_perf,
const std::vector<double>& z_perf,
const std::vector<double>& surf_dens,
const double gravity)
Opm::WellDensitySegmented::computeConnectionDensities(const Wells& wells,
const WellStateFullyImplicitBlackoil& wstate,
const PhaseUsage& phase_usage,
const std::vector<double>& b_perf,
const std::vector<double>& rsmax_perf,
const std::vector<double>& rvmax_perf,
const std::vector<double>& surf_dens)
{
// Verify that we have consistent input.
const int np = wells.number_of_phases;
@@ -53,9 +54,6 @@ Opm::WellDensitySegmented::computeConnectionPressureDelta(const Wells& wells,
if (nperf*np != int(b_perf.size())) {
OPM_THROW(std::logic_error, "Inconsistent input: wells vs. b_perf.");
}
if (nperf != int(z_perf.size())) {
OPM_THROW(std::logic_error, "Inconsistent input: wells vs. z_perf.");
}
if ((!rsmax_perf.empty()) || (!rvmax_perf.empty())) {
// Need both oil and gas phases.
if (!phase_usage.phase_used[BlackoilPhases::Liquid]) {
@@ -66,14 +64,6 @@ Opm::WellDensitySegmented::computeConnectionPressureDelta(const Wells& wells,
}
}
// Algorithm:
// We'll assume the perforations are given in order from top to
// bottom for each well. By top and bottom we do not necessarily
// mean in a geometric sense (depth), but in a topological sense:
// the 'top' perforation is nearest to the surface topologically.
// Our goal is to compute a pressure delta for each perforation.
// 1. Compute the flow (in surface volume units for each
// component) exiting up the wellbore from each perforation,
// taking into account flow from lower in the well, and
@@ -145,7 +135,37 @@ Opm::WellDensitySegmented::computeConnectionPressureDelta(const Wells& wells,
}
}
// 3. Compute pressure differences between perforations.
return dens;
}
std::vector<double>
Opm::WellDensitySegmented::computeConnectionPressureDelta(const Wells& wells,
const std::vector<double>& z_perf,
const std::vector<double>& dens_perf,
const double gravity) {
const int nw = wells.number_of_wells;
const int nperf = wells.well_connpos[nw];
if (nperf != int(z_perf.size())) {
OPM_THROW(std::logic_error, "Inconsistent input: wells vs. z_perf.");
}
if (nperf != int(dens_perf.size())) {
OPM_THROW(std::logic_error, "Inconsistent input: wells vs. dens_perf.");
}
// Algorithm:
// We'll assume the perforations are given in order from top to
// bottom for each well. By top and bottom we do not necessarily
// mean in a geometric sense (depth), but in a topological sense:
// the 'top' perforation is nearest to the surface topologically.
// Our goal is to compute a pressure delta for each perforation.
// 1. Compute pressure differences between perforations.
// dp_perf will contain the pressure difference between a
// perforation and the one above it, except for the first
// perforation for each well, for which it will be the
@@ -155,11 +175,11 @@ Opm::WellDensitySegmented::computeConnectionPressureDelta(const Wells& wells,
for (int perf = wells.well_connpos[w]; perf < wells.well_connpos[w+1]; ++perf) {
const double z_above = perf == wells.well_connpos[w] ? wells.depth_ref[w] : z_perf[perf - 1];
const double dz = z_perf[perf] - z_above;
dp_perf[perf] = dz * dens[perf] * gravity;
dp_perf[perf] = dz * dens_perf[perf] * gravity;
}
}
// 4. Compute pressure differences to the reference point (bhp) by
// 2. Compute pressure differences to the reference point (bhp) by
// accumulating the already computed adjacent pressure
// differences, storing the result in dp_perf.
// This accumulation must be done per well.

23
opm/autodiff/WellDensitySegmented.hpp
View File

@@ -39,7 +39,7 @@ namespace Opm
class WellDensitySegmented
{
public:
/// Compute pressure deltas.
/// Compute well segment densities
/// Notation: N = number of perforations, P = number of phases.
/// \param[in] wells struct with static well info
/// \param[in] wstate dynamic well solution information, only perfRates() is used
@@ -47,17 +47,24 @@ namespace Opm
/// \param[in] b_perf inverse ('little b') formation volume factor, size NP, P values per perforation
/// \param[in] rsmax_perf saturation point for rs (gas in oil) at each perforation, size N
/// \param[in] rvmax_perf saturation point for rv (oil in gas) at each perforation, size N
/// \param[in] z_perf depth values for each perforation, size N
/// \param[in] surf_dens surface densities for active components, size P
static std::vector<double> computeConnectionDensities(const Wells& wells,
const WellStateFullyImplicitBlackoil& wstate,
const PhaseUsage& phase_usage,
const std::vector<double>& b_perf,
const std::vector<double>& rsmax_perf,
const std::vector<double>& rvmax_perf,
const std::vector<double>& surf_dens);
/// Compute pressure deltas.
/// Notation: N = number of perforations, P = number of phases.
/// \param[in] wells struct with static well info
/// \param[in] z_perf depth values for each perforation, size N
/// \param[in] dens_perf densities for each perforation, size N (typically computed using computeConnectionDensities)
/// \param[in] gravity gravity acceleration constant
static std::vector<double> computeConnectionPressureDelta(const Wells& wells,
const WellStateFullyImplicitBlackoil& wstate,
const PhaseUsage& phase_usage,
const std::vector<double>& b_perf,
const std::vector<double>& rsmax_perf,
const std::vector<double>& rvmax_perf,
const std::vector<double>& z_perf,
const std::vector<double>& surf_dens,
const std::vector<double>& dens_perf,
const double gravity);
};

2306
tests/VFPPROD1 Normal file
View File

File diff suppressed because it is too large Load Diff

2295
tests/VFPPROD2 Normal file
View File

File diff suppressed because it is too large Load Diff

72
tests/test_autodiffhelpers.cpp
View File

@@ -144,3 +144,75 @@ BOOST_AUTO_TEST_CASE(vertcatCollapseJacsTest)
BOOST_CHECK((x.value() == expected_val).all());
BOOST_CHECK(x.derivative()[0] == expected_jac);
}
BOOST_AUTO_TEST_CASE(supersetTest)
{
typedef AutoDiffBlock<double> ADB;
{
ADB subset = ADB::constant(ADB::V::Ones(3));
const int full_size = 32;
std::vector<int> indices = {1, 3, 5};
AutoDiffBlock<double> n_vals = superset(subset, indices, full_size);
BOOST_CHECK_EQUAL(n_vals.size(), full_size);
for (int i=0; i<n_vals.size(); ++i) {
if (find(indices.begin(), indices.end(), i) != indices.end()) {
BOOST_CHECK_EQUAL(n_vals.value()[i], 1);
}
else {
BOOST_CHECK_EQUAL(n_vals.value()[i], 0);
}
}
}
}
BOOST_AUTO_TEST_CASE(supersetEmptyTest)
{
typedef AutoDiffBlock<double> ADB;
{
ADB subset = ADB::constant(ADB::V::Ones(0));
const int full_size = 32;
std::vector<int> indices = {};
AutoDiffBlock<double> n_vals = superset(subset, indices, full_size);
BOOST_CHECK_EQUAL(n_vals.size(), full_size);
for (int i=0; i<n_vals.size(); ++i) {
BOOST_CHECK_EQUAL(n_vals.value()[i], 0);
}
}
}
BOOST_AUTO_TEST_CASE(subsetTest)
{
typedef AutoDiffBlock<double> ADB;
{
ADB superset = ADB::constant(ADB::V::Ones(32));
std::vector<int> indices = {1, 3, 5};
AutoDiffBlock<double> n_vals = subset(superset, indices);
BOOST_CHECK_EQUAL(n_vals.size(), 3);
for (int i=0; i<n_vals.size(); ++i) {
if (find(indices.begin(), indices.end(), i) != indices.end()) {
BOOST_CHECK_EQUAL(n_vals.value()[i], 1);
}
}
}
}
BOOST_AUTO_TEST_CASE(subsetEmptyTest)
{
typedef AutoDiffBlock<double> ADB;
{
ADB superset = ADB::constant(ADB::V::Ones(32));
std::vector<int> indices = {};
AutoDiffBlock<double> n_vals = subset(superset, indices);
BOOST_CHECK_EQUAL(n_vals.size(), 0);
}
}

1162
tests/test_vfpproperties.cpp Normal file
View File

File diff suppressed because one or more lines are too long

10
tests/test_welldensitysegmented.cpp
View File

@@ -90,7 +90,15 @@ BOOST_AUTO_TEST_CASE(TestPressureDeltas)
const std::vector<double> z_perf = { 10, 30, 50, 70, 90, 10, 30, 50, 70, 90 };
const std::vector<double> surf_dens = { 1000.0, 800.0, 10.0 };
const double gravity = Opm::unit::gravity;
const std::vector<double> dp = WellDensitySegmented::computeConnectionPressureDelta(*wells, wellstate, pu, b_perf, rsmax_perf, rvmax_perf, z_perf, surf_dens, gravity);
std::vector<double> cd =
WellDensitySegmented::computeConnectionDensities(
*wells, wellstate, pu,
b_perf, rsmax_perf, rvmax_perf, surf_dens);
std::vector<double> dp =
WellDensitySegmented::computeConnectionPressureDelta(
*wells, z_perf, cd, gravity);
const std::vector<double> answer = { 20e3*gravity, 62e3*gravity, 106e3*gravity, 152e3*gravity, 200e3*gravity,
20e3*gravity, 62e3*gravity, 106e3*gravity, 152e3*gravity, 200e3*gravity };
BOOST_REQUIRE_EQUAL(dp.size(), answer.size());