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 #5353 from akva2/vfp_template_scalar

Browse Source 
VFP classes: template Scalar type
This commit is contained in:
Bård Skaflestad
2024-05-21 17:52:10 +02:00
committed by GitHub
parent 15001cbf7e d5d16eaee4
commit 9b1797e7b7
16 changed files with 558 additions and 498 deletions
Download Patch File Download Diff File Expand all files Collapse all files

18
examples/printvfp.cpp
View File

@@ -46,7 +46,7 @@ struct Setup
std::shared_ptr<Python> python;
std::unique_ptr<const Schedule> schedule;
std::unique_ptr<SummaryState> summary_state;
std::unique_ptr<VFPProperties> vfp_properties;
std::unique_ptr<VFPProperties<double>> vfp_properties;
Setup(const std::string& file)
{
@@ -63,7 +63,9 @@ struct Setup
const int step = 0;
const auto& sched_state = schedule->operator[](step);
WellState<double> well_state(phaseUsage(runspec.phases()));
vfp_properties = std::make_unique<VFPProperties>(sched_state.vfpinj(), sched_state.vfpprod(), well_state);
vfp_properties = std::make_unique<VFPProperties<double>>(sched_state.vfpinj(),
sched_state.vfpprod(),
well_state);
};
};
@@ -81,13 +83,13 @@ double computeBhp(const VFPProdTable& table,
// First, find the values to interpolate between.
// Assuming positive flo here!
assert(flo > 0.0);
auto flo_i = detail::findInterpData(flo, table.getFloAxis());
auto thp_i = detail::findInterpData(thp, table.getTHPAxis()); // assume constant
auto wfr_i = detail::findInterpData(wfr, table.getWFRAxis());
auto gfr_i = detail::findInterpData(gfr, table.getGFRAxis());
auto alq_i = detail::findInterpData(alq, table.getALQAxis()); //assume constant
auto flo_i = VFPHelpers<double>::findInterpData(flo, table.getFloAxis());
auto thp_i = VFPHelpers<double>::findInterpData(thp, table.getTHPAxis()); // assume constant
auto wfr_i = VFPHelpers<double>::findInterpData(wfr, table.getWFRAxis());
auto gfr_i = VFPHelpers<double>::findInterpData(gfr, table.getGFRAxis());
auto alq_i = VFPHelpers<double>::findInterpData(alq, table.getALQAxis()); //assume constant
return detail::interpolate(table, flo_i, thp_i, wfr_i, gfr_i, alq_i).value;
return VFPHelpers<double>::interpolate(table, flo_i, thp_i, wfr_i, gfr_i, alq_i).value;
}

4
opm/simulators/wells/BlackoilWellModelGeneric.hpp
View File

@@ -63,7 +63,7 @@ namespace Opm {
class Schedule;
struct SimulatorUpdate;
class SummaryConfig;
class VFPProperties;
template<class Scalar> class VFPProperties;
template<class Scalar> class WellInterfaceGeneric;
template<class Scalar> class WellState;
} // namespace Opm
@@ -557,7 +557,7 @@ protected:
mutable std::unordered_set<std::string> closed_this_step_;
GuideRate guideRate_;
std::unique_ptr<VFPProperties> vfp_properties_{};
std::unique_ptr<VFPProperties<Scalar>> vfp_properties_{};
std::map<std::string, Scalar> node_pressures_; // Storing network pressures for output.
// previous injection multiplier, it is used in the injection multiplier calculation for WINJMULT keyword

2
opm/simulators/wells/BlackoilWellModel_impl.hpp
View File

@@ -275,7 +275,7 @@ namespace Opm {
{
const auto& sched_state = this->schedule()[timeStepIdx];
this->vfp_properties_ = std::make_unique<VFPProperties>
this->vfp_properties_ = std::make_unique<VFPProperties<Scalar>>
(sched_state.vfpinj(), sched_state.vfpprod(), this->wellState());
}
}

362
opm/simulators/wells/VFPHelpers.cpp
View File

@@ -38,14 +38,15 @@ namespace {
* Helper function that finds x for a given value of y for a line
* *NOTE ORDER OF ARGUMENTS*
*/
double findX(const double x0,
const double x1,
const double y0,
const double y1,
const double y)
template<class Scalar>
Scalar findX(const Scalar x0,
const Scalar x1,
const Scalar y0,
const Scalar y1,
const Scalar y)
{
const double dx = x1 - x0;
const double dy = y1 - y0;
const Scalar dx = x1 - x0;
const Scalar dy = y1 - y0;
/**
* y = y0 + (dy / dx) * (x - x0)
@@ -54,7 +55,7 @@ double findX(const double x0,
* If dy is zero, use x1 as the value.
*/
double x = 0.0;
Scalar x = 0.0;
if (dy != 0.0) {
x = x0 + (y-y0) * (dx/dy);
@@ -77,17 +78,18 @@ static T chopNegativeValues(const T& value) {
}
namespace Opm {
namespace detail {
InterpData findInterpData(const double value_in, const std::vector<double>& values)
template<class Scalar>
detail::InterpData<Scalar> VFPHelpers<Scalar>::findInterpData(const Scalar value_in,
const std::vector<double>& values)
{
InterpData retval;
detail::InterpData<Scalar> retval;
const int nvalues = values.size();
// chopping the value to be zero, which means we do not
// extrapolate the table towards nagative ranges
const double value = value_in < 0.? 0. : value_in;
const Scalar value = value_in < 0.? 0. : value_in;
//If we only have one value in our vector, return that
if (values.size() == 1) {
@@ -119,8 +121,8 @@ InterpData findInterpData(const double value_in, const std::vector<double>& valu
}
}
const double start = values[retval.ind_[0]];
const double end = values[retval.ind_[1]];
const Scalar start = values[retval.ind_[0]];
const Scalar end = values[retval.ind_[1]];
//Find interpolation ratio
if (end > start) {
@@ -138,50 +140,17 @@ InterpData findInterpData(const double value_in, const std::vector<double>& valu
return retval;
}
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;
}
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;
}
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;
}
VFPEvaluation interpolate(const VFPProdTable& table,
const InterpData& flo_i,
const InterpData& thp_i,
const InterpData& wfr_i,
const InterpData& gfr_i,
const InterpData& alq_i)
template<class Scalar>
detail::VFPEvaluation<Scalar> VFPHelpers<Scalar>::
interpolate(const VFPProdTable& table,
const detail::InterpData<Scalar>& flo_i,
const detail::InterpData<Scalar>& thp_i,
const detail::InterpData<Scalar>& wfr_i,
const detail::InterpData<Scalar>& gfr_i,
const detail::InterpData<Scalar>& alq_i)
{
//Values and derivatives in a 5D hypercube
VFPEvaluation nn[2][2][2][2][2];
detail::VFPEvaluation<Scalar> 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
@@ -231,7 +200,7 @@ VFPEvaluation interpolate(const VFPProdTable& table,
}
}
double t1, t2; //interpolation variables, so that t1 = (1-t) and t2 = t.
Scalar 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
@@ -280,13 +249,14 @@ VFPEvaluation interpolate(const VFPProdTable& table,
return nn[0][0][0][0][0];
}
VFPEvaluation interpolate(const VFPInjTable& table,
const InterpData& flo_i,
const InterpData& thp_i)
template<class Scalar>
detail::VFPEvaluation<Scalar> VFPHelpers<Scalar>::
interpolate(const VFPInjTable& table,
const detail::InterpData<Scalar>& flo_i,
const detail::InterpData<Scalar>& thp_i)
{
//Values and derivatives in a 2D plane
VFPEvaluation nn[2][2];
detail::VFPEvaluation<Scalar> 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.
@@ -318,7 +288,7 @@ VFPEvaluation interpolate(const VFPInjTable& table,
nn[i][1].dflo = nn[i][0].dflo;
}
double t1, t2; //interpolation variables, so that t1 = (1-t) and t2 = t.
Scalar t1, t2; //interpolation variables, so that t1 = (1-t) and t2 = t.
// Go from 2D to 1D
t2 = flo_i.factor_;
@@ -334,20 +304,22 @@ VFPEvaluation interpolate(const VFPInjTable& table,
return nn[0][0];
}
VFPEvaluation bhp(const VFPProdTable& table,
const double aqua,
const double liquid,
const double vapour,
const double thp,
const double alq,
const double explicit_wfr,
const double explicit_gfr,
const bool use_vfpexplicit)
template<class Scalar>
detail::VFPEvaluation<Scalar> VFPHelpers<Scalar>::
bhp(const VFPProdTable& table,
const Scalar aqua,
const Scalar liquid,
const Scalar vapour,
const Scalar thp,
const Scalar alq,
const Scalar explicit_wfr,
const Scalar explicit_gfr,
const bool use_vfpexplicit)
{
//Find interpolation variables
double flo = detail::getFlo(table, aqua, liquid, vapour);
double wfr = detail::getWFR(table, aqua, liquid, vapour);
double gfr = detail::getGFR(table, aqua, liquid, vapour);
Scalar flo = detail::getFlo(table, aqua, liquid, vapour);
Scalar wfr = detail::getWFR(table, aqua, liquid, vapour);
Scalar gfr = detail::getGFR(table, aqua, liquid, vapour);
if (use_vfpexplicit || -flo < table.getFloAxis().front()) {
wfr = explicit_wfr;
gfr = explicit_gfr;
@@ -355,43 +327,47 @@ VFPEvaluation bhp(const VFPProdTable& table,
//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());
auto flo_i = findInterpData(-flo, table.getFloAxis());
auto thp_i = findInterpData( thp, table.getTHPAxis());
auto wfr_i = findInterpData( wfr, table.getWFRAxis());
auto gfr_i = findInterpData( gfr, table.getGFRAxis());
auto alq_i = findInterpData( alq, table.getALQAxis());
detail::VFPEvaluation retval = detail::interpolate(table, flo_i, thp_i, wfr_i, gfr_i, alq_i);
detail::VFPEvaluation retval = interpolate(table, flo_i, thp_i, wfr_i, gfr_i, alq_i);
return retval;
}
VFPEvaluation bhp(const VFPInjTable& table,
const double aqua,
const double liquid,
const double vapour,
const double thp)
template<class Scalar>
detail::VFPEvaluation<Scalar> VFPHelpers<Scalar>::
bhp(const VFPInjTable& table,
const Scalar aqua,
const Scalar liquid,
const Scalar vapour,
const Scalar thp)
{
//Find interpolation variables
double flo = detail::getFlo(table, aqua, liquid, vapour);
Scalar flo = detail::getFlo(table, aqua, liquid, vapour);
//First, find the values to interpolate between
auto flo_i = detail::findInterpData(flo, table.getFloAxis());
auto thp_i = detail::findInterpData(thp, table.getTHPAxis());
auto flo_i = findInterpData(flo, table.getFloAxis());
auto thp_i = findInterpData(thp, table.getTHPAxis());
//Then perform the interpolation itself
detail::VFPEvaluation retval = detail::interpolate(table, flo_i, thp_i);
detail::VFPEvaluation retval = interpolate(table, flo_i, thp_i);
return retval;
}
double findTHP(const std::vector<double>& bhp_array,
const std::vector<double>& thp_array,
double bhp)
template<class Scalar>
Scalar VFPHelpers<Scalar>::
findTHP(const std::vector<Scalar>& bhp_array,
const std::vector<double>& thp_array,
Scalar bhp)
{
int nthp = thp_array.size();
double thp = -1e100;
Scalar thp = -1e100;
//Check that our thp axis is sorted
assert(std::is_sorted(thp_array.begin(), thp_array.end()));
@@ -406,19 +382,19 @@ double findTHP(const std::vector<double>& bhp_array,
//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];
const Scalar& x0 = thp_array[0];
const Scalar& x1 = thp_array[1];
const Scalar& y0 = bhp_array[0];
const Scalar& y1 = bhp_array[1];
thp = 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];
const Scalar& x0 = thp_array[nthp-2];
const Scalar& x1 = thp_array[nthp-1];
const Scalar& y0 = bhp_array[nthp-2];
const Scalar& y1 = bhp_array[nthp-1];
thp = findX(x0, x1, y0, y1, bhp);
}
//Target bhp within table ranges, interpolate
@@ -432,8 +408,8 @@ double findTHP(const std::vector<double>& bhp_array,
int i=0;
bool found = false;
for (; i<nthp-1; ++i) {
const double& y0 = bhp_array[i ];
const double& y1 = bhp_array[i+1];
const Scalar& y0 = bhp_array[i ];
const Scalar& y1 = bhp_array[i+1];
if (y0 < bhp && bhp <= y1) {
found = true;
@@ -444,10 +420,10 @@ double findTHP(const std::vector<double>& bhp_array,
assert(found == true);
static_cast<void>(found); //Silence compiler warning
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];
const Scalar& x0 = thp_array[i ];
const Scalar& x1 = thp_array[i+1];
const Scalar& y0 = bhp_array[i ];
const Scalar& y1 = bhp_array[i+1];
thp = findX(x0, x1, y0, y1, bhp);
}
}
@@ -459,8 +435,8 @@ double findTHP(const std::vector<double>& bhp_array,
int i=0;
bool found = false;
for (; i<nthp-1; ++i) {
const double& y0 = bhp_array[i ];
const double& y1 = bhp_array[i+1];
const Scalar& y0 = bhp_array[i ];
const Scalar& y1 = bhp_array[i+1];
if (y0 < bhp && bhp <= y1) {
found = true;
@@ -468,27 +444,27 @@ double findTHP(const std::vector<double>& bhp_array,
}
}
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];
const Scalar& x0 = thp_array[i ];
const Scalar& x1 = thp_array[i+1];
const Scalar& y0 = bhp_array[i ];
const Scalar& y1 = bhp_array[i+1];
thp = 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];
const Scalar& x0 = thp_array[0];
const Scalar& x1 = thp_array[1];
const Scalar& y0 = bhp_array[0];
const Scalar& y1 = bhp_array[1];
thp = 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];
const Scalar& x0 = thp_array[nthp-2];
const Scalar& x1 = thp_array[nthp-1];
const Scalar& y0 = bhp_array[nthp-2];
const Scalar& y1 = bhp_array[nthp-1];
thp = findX(x0, x1, y0, y1, bhp);
}
else {
@@ -499,86 +475,88 @@ double findTHP(const std::vector<double>& bhp_array,
return thp;
}
std::pair<double, double>
template<class Scalar>
std::pair<Scalar, Scalar> VFPHelpers<Scalar>::
getMinimumBHPCoordinate(const VFPProdTable& table,
const double thp,
const double wfr,
const double gfr,
const double alq)
{
const Scalar thp,
const Scalar wfr,
const Scalar gfr,
const Scalar alq)
{
// Given fixed thp, wfr, gfr and alq, this function finds the minimum bhp and returns
// the corresponding pair (-flo_at_bhp_min, bhp_min). No assumption is taken on the
// shape of the function bhp(flo), so all points in the flo-axis is checked.
double flo_at_bhp_min = 0.0; // start by checking flo=0
auto flo_i = detail::findInterpData(flo_at_bhp_min, 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());
// the corresponding pair (-flo_at_bhp_min, bhp_min). No assumption is taken on the
// shape of the function bhp(flo), so all points in the flo-axis is checked.
Scalar flo_at_bhp_min = 0.0; // start by checking flo=0
auto flo_i = findInterpData(flo_at_bhp_min, table.getFloAxis());
auto thp_i = findInterpData( thp, table.getTHPAxis());
auto wfr_i = findInterpData( wfr, table.getWFRAxis());
auto gfr_i = findInterpData( gfr, table.getGFRAxis());
auto alq_i = findInterpData( alq, table.getALQAxis());
detail::VFPEvaluation bhp_i = detail::interpolate(table, flo_i, thp_i, wfr_i, gfr_i, alq_i);
double bhp_min = bhp_i.value;
detail::VFPEvaluation bhp_i = interpolate(table, flo_i, thp_i, wfr_i, gfr_i, alq_i);
Scalar bhp_min = bhp_i.value;
const std::vector<double>& flos = table.getFloAxis();
for (size_t i = 0; i < flos.size(); ++i) {
flo_i = detail::findInterpData(flos[i], flos);
bhp_i = detail::interpolate(table, flo_i, thp_i, wfr_i, gfr_i, alq_i);
flo_i = findInterpData(flos[i], flos);
bhp_i = interpolate(table, flo_i, thp_i, wfr_i, gfr_i, alq_i);
if (bhp_i.value < bhp_min){
bhp_min = bhp_i.value;
flo_at_bhp_min = flos[i];
}
}
// return negative flo
return std::make_pair(-flo_at_bhp_min, bhp_min);
}
return std::make_pair(-flo_at_bhp_min, bhp_min);
}
std::optional<std::pair<double, double>>
template<class Scalar>
std::optional<std::pair<Scalar, Scalar>> VFPHelpers<Scalar>::
intersectWithIPR(const VFPProdTable& table,
const double thp,
const double wfr,
const double gfr,
const double alq,
const double ipr_a,
const double ipr_b,
const std::function<double(const double)>& adjust_bhp)
{
// Given fixed thp, wfr, gfr and alq, this function finds a stable (-flo, bhp)-intersection
// between the ipr-line and bhp(flo) if such an intersection exists. For multiple stable
const Scalar thp,
const Scalar wfr,
const Scalar gfr,
const Scalar alq,
const Scalar ipr_a,
const Scalar ipr_b,
const std::function<Scalar(const Scalar)>& adjust_bhp)
{
// Given fixed thp, wfr, gfr and alq, this function finds a stable (-flo, bhp)-intersection
// between the ipr-line and bhp(flo) if such an intersection exists. For multiple stable
// intersections, the one corresponding the largest flo is returned.
// The adjust_bhp-function is used to adjust the vfp-table bhp-values to actual bhp-values due
// vfp/well ref-depth differences and/or WVFPDP-related pressure adjustments.
// vfp/well ref-depth differences and/or WVFPDP-related pressure adjustments.
// NOTE: ipr-line is q=b*bhp - a!
// ipr is given for negative flo, so
// flo = -b*bhp + a, i.e., bhp = -(flo-a)/b
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());
// flo = -b*bhp + a, i.e., bhp = -(flo-a)/b
auto thp_i = findInterpData( thp, table.getTHPAxis());
auto wfr_i = findInterpData( wfr, table.getWFRAxis());
auto gfr_i = findInterpData( gfr, table.getGFRAxis());
auto alq_i = findInterpData( alq, table.getALQAxis());
if (ipr_b == 0.0) {
// this shouldn't happen, but deal with it to be safe
auto flo_i = detail::findInterpData(ipr_a, table.getFloAxis());
detail::VFPEvaluation bhp_i = detail::interpolate(table, flo_i, thp_i, wfr_i, gfr_i, alq_i);
auto flo_i = findInterpData(ipr_a, table.getFloAxis());
detail::VFPEvaluation bhp_i = interpolate(table, flo_i, thp_i, wfr_i, gfr_i, alq_i);
return std::make_pair(-ipr_a, adjust_bhp(bhp_i.value));
}
// find largest flo (flo_x) for which y = bhp(flo) + (flo-a)/b = 0 and dy/dflo > 0
double flo_x = -1.0;
double flo0, flo1;
double y0, y1;
Scalar flo_x = -1.0;
Scalar flo0, flo1;
Scalar y0, y1;
flo0 = 0.0; // start by checking flo=0
auto flo_i = detail::findInterpData(flo0, table.getFloAxis());
detail::VFPEvaluation bhp_i = detail::interpolate(table, flo_i, thp_i, wfr_i, gfr_i, alq_i);
auto flo_i = findInterpData(flo0, table.getFloAxis());
detail::VFPEvaluation bhp_i = interpolate(table, flo_i, thp_i, wfr_i, gfr_i, alq_i);
y0 = adjust_bhp(bhp_i.value) - ipr_a/ipr_b; // +0.0/ipr_b
std::vector<double> flos = table.getFloAxis();
const std::vector<double>& flos = table.getFloAxis();
for (size_t i = 0; i < flos.size(); ++i) {
flo1 = flos[i];
flo_i = detail::findInterpData(flo1, table.getFloAxis());
bhp_i = detail::interpolate(table, flo_i, thp_i, wfr_i, gfr_i, alq_i);
flo_i = findInterpData(flo1, flos);
bhp_i = interpolate(table, flo_i, thp_i, wfr_i, gfr_i, alq_i);
y1 = adjust_bhp(bhp_i.value) + (flo1 - ipr_a)/ipr_b;
if (y0 < 0 && y1 >= 0){
// crossing with positive slope
double w = -y0/(y1-y0);
Scalar w = -y0/(y1-y0);
w = std::clamp(w, 0.0, 1.0); // just to be safe (if y0~y1~0)
flo_x = flo0 + w*(flo1 - flo0);
}
@@ -591,7 +569,46 @@ intersectWithIPR(const VFPProdTable& table,
} else {
return std::nullopt;
}
}
}
namespace detail {
template<class Scalar>
VFPEvaluation<Scalar> operator+(VFPEvaluation<Scalar> lhs, const VFPEvaluation<Scalar>& 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;
}
template<class Scalar>
VFPEvaluation<Scalar> operator-(VFPEvaluation<Scalar> lhs, const VFPEvaluation<Scalar>& 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;
}
template<class Scalar>
VFPEvaluation<Scalar> operator*(Scalar lhs, const VFPEvaluation<Scalar>& rhs)
{
VFPEvaluation<Scalar> 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;
}
template <typename T>
T getFlo(const VFPProdTable& table,
@@ -753,4 +770,7 @@ INSTANCE(DenseAd::Evaluation<double, 9, 0u>)
INSTANCE(DenseAd::Evaluation<double, 10, 0u>)
} // namespace detail
template class VFPHelpers<double>;
} // namespace Opm

218
opm/simulators/wells/VFPHelpers.hpp
View File

@@ -21,7 +21,6 @@
#ifndef OPM_AUTODIFF_VFPHELPERS_HPP_
#define OPM_AUTODIFF_VFPHELPERS_HPP_
#include <array>
#include <functional>
#include <map>
#include <vector>
@@ -37,6 +36,40 @@ class VFPProdTable;
namespace detail {
/**
* An "ADB-like" structure with a single value and a set of derivatives
*/
template<class Scalar>
struct VFPEvaluation
{
VFPEvaluation() : value(0.0), dthp(0.0), dwfr(0.0), dgfr(0.0), dalq(0.0), dflo(0.0) {};
Scalar value;
Scalar dthp;
Scalar dwfr;
Scalar dgfr;
Scalar dalq;
Scalar dflo;
};
template<class Scalar>
VFPEvaluation<Scalar> operator+(VFPEvaluation<Scalar> lhs, const VFPEvaluation<Scalar>& rhs);
template<class Scalar>
VFPEvaluation<Scalar> operator-(VFPEvaluation<Scalar> lhs, const VFPEvaluation<Scalar>& rhs);
template<class Scalar>
VFPEvaluation<Scalar> operator*(Scalar lhs, const VFPEvaluation<Scalar>& rhs);
/**
* Helper struct for linear interpolation
*/
template<class Scalar>
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]
Scalar 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
Scalar factor_; // Interpolation factor
};
/**
* Computes the flo parameter according to the flo_type_
* for production tables
@@ -79,76 +112,6 @@ T getGFR(const VFPProdTable& table,
const T& liquid,
const T& vapour);
/**
* 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_in Value to find in values
* @param values Sorted list of values to search for value in.
* @return Data required to find the interpolated value
*/
InterpData findInterpData(const double value_in, const std::vector<double>& values);
/**
* 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;
};
VFPEvaluation operator+(VFPEvaluation lhs, const VFPEvaluation& rhs);
VFPEvaluation operator-(VFPEvaluation lhs, const VFPEvaluation& rhs);
VFPEvaluation operator*(double lhs, const VFPEvaluation& rhs);
/**
* Helper function which interpolates data using the indices etc. given in the inputs.
*/
VFPEvaluation interpolate(const VFPProdTable& table,
const InterpData& flo_i,
const InterpData& thp_i,
const InterpData& wfr_i,
const InterpData& gfr_i,
const InterpData& alq_i);
/**
* This basically models interpolate(VFPProdTable::array_type, ...)
* which performs 5D interpolation, but here for the 2D case only
*/
VFPEvaluation interpolate(const VFPInjTable& table,
const InterpData& flo_i,
const InterpData& thp_i);
VFPEvaluation bhp(const VFPProdTable& table,
const double aqua,
const double liquid,
const double vapour,
const double thp,
const double alq,
const double explicit_wfr,
const double explicit_gfr,
const bool use_vfpexplicit);
VFPEvaluation bhp(const VFPInjTable& table,
const double aqua,
const double liquid,
const double vapour,
const double thp);
/**
* Returns the table from the map if found, or throws an exception
*/
@@ -164,54 +127,95 @@ bool hasTable(const std::map<int, std::reference_wrapper<const T>>& tables, int
return (entry != tables.end() );
}
/**
* Returns the type variable for FLO/GFR/WFR for production tables
*/
template <typename TYPE, typename TABLE>
TYPE getType(const TABLE& table);
}
/**
* 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.
*/
double findTHP(const std::vector<double>& bhp_array,
const std::vector<double>& thp_array,
double bhp);
template<class Scalar>
class VFPHelpers {
public:
/**
* Helper function to find indices etc. for linear interpolation and extrapolation
* @param value_in Value to find in values
* @param values Sorted list of values to search for value in.
* @return Data required to find the interpolated value
*/
static detail::InterpData<Scalar> findInterpData(const Scalar value_in,
const std::vector<double>& values);
/**
* Get (flo, bhp) at minimum bhp for given thp,wfr,gfr,alq
*/
std::pair<double, double>
getMinimumBHPCoordinate(const VFPProdTable& table,
const double thp,
const double wfr,
const double gfr,
const double alq);
/**
* Helper function which interpolates data using the indices etc. given in the inputs.
*/
static detail::VFPEvaluation<Scalar> interpolate(const VFPProdTable& table,
const detail::InterpData<Scalar>& flo_i,
const detail::InterpData<Scalar>& thp_i,
const detail::InterpData<Scalar>& wfr_i,
const detail::InterpData<Scalar>& gfr_i,
const detail::InterpData<Scalar>& alq_i);
/**
* Get (flo, bhp) at largest occuring stable vfp/ipr-intersection
* if it exists
*/
std::optional<std::pair<double, double>>
intersectWithIPR(const VFPProdTable& table,
const double thp,
const double wfr,
const double gfr,
const double alq,
const double ipr_a,
const double ipr_b,
const std::function<double(const double)>& adjust_bhp);
/**
* This basically models interpolate(VFPProdTable::array_type, ...)
* which performs 5D interpolation, but here for the 2D case only
*/
static detail::VFPEvaluation<Scalar> interpolate(const VFPInjTable& table,
const detail::InterpData<Scalar>& flo_i,
const detail::InterpData<Scalar>& thp_i);
} // namespace detail
static detail::VFPEvaluation<Scalar> bhp(const VFPProdTable& table,
const Scalar aqua,
const Scalar liquid,
const Scalar vapour,
const Scalar thp,
const Scalar alq,
const Scalar explicit_wfr,
const Scalar explicit_gfr,
const bool use_vfpexplicit);
static detail::VFPEvaluation<Scalar> bhp(const VFPInjTable& table,
const Scalar aqua,
const Scalar liquid,
const Scalar vapour,
const Scalar thp);
/**
* 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.
*/
static Scalar findTHP(const std::vector<Scalar>& bhp_array,
const std::vector<double>& thp_array,
Scalar bhp);
/**
* Get (flo, bhp) at minimum bhp for given thp,wfr,gfr,alq
*/
static std::pair<Scalar, Scalar>
getMinimumBHPCoordinate(const VFPProdTable& table,
const Scalar thp,
const Scalar wfr,
const Scalar gfr,
const Scalar alq);
/**
* Get (flo, bhp) at largest occuring stable vfp/ipr-intersection
* if it exists
*/
static std::optional<std::pair<Scalar, Scalar>>
intersectWithIPR(const VFPProdTable& table,
const Scalar thp,
const Scalar wfr,
const Scalar gfr,
const Scalar alq,
const Scalar ipr_a,
const Scalar ipr_b,
const std::function<Scalar(const Scalar)>& adjust_bhp);
};
} // namespace
#endif /* OPM_AUTODIFF_VFPHELPERS_HPP_ */

98
opm/simulators/wells/VFPInjProperties.cpp
View File

@@ -17,49 +17,51 @@
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#include "config.h"
#include <opm/simulators/wells/VFPInjProperties.hpp>
#include <opm/input/eclipse/Schedule/VFPInjTable.hpp>
#include <opm/material/densead/Math.hpp>
#include <opm/material/densead/Evaluation.hpp>
#include <opm/input/eclipse/Schedule/VFPInjTable.hpp>
#include <opm/simulators/wells/VFPHelpers.hpp>
#include <limits>
namespace Opm {
double VFPInjProperties::bhp(int table_id,
const double& aqua,
const double& liquid,
const double& vapour,
const double& thp_arg) const {
template<class Scalar>
Scalar VFPInjProperties<Scalar>::
bhp(const int table_id,
const Scalar aqua,
const Scalar liquid,
const Scalar vapour,
const Scalar thp_arg) const
{
const VFPInjTable& table = detail::getTable(m_tables, table_id);
detail::VFPEvaluation retval = detail::bhp(table, aqua, liquid, vapour, thp_arg);
detail::VFPEvaluation retval = VFPHelpers<Scalar>::bhp(table, aqua, liquid, vapour, thp_arg);
return retval.value;
}
double VFPInjProperties::thp(int table_id,
const double& aqua,
const double& liquid,
const double& vapour,
const double& bhp_arg) const {
template<class Scalar>
Scalar VFPInjProperties<Scalar>::
thp(const int table_id,
const Scalar aqua,
const Scalar liquid,
const Scalar vapour,
const Scalar bhp_arg) const
{
const VFPInjTable& table = detail::getTable(m_tables, table_id);
//Find interpolation variables
const double flo = detail::getFlo(table, aqua, liquid, vapour);
if (std::abs(flo) < std::numeric_limits<double>::epsilon()) {
const Scalar flo = detail::getFlo(table, aqua, liquid, vapour);
if (std::abs(flo) < std::numeric_limits<Scalar>::epsilon()) {
return 0.;
}
const std::vector<double> thp_array = table.getTHPAxis();
const auto thp_array = table.getTHPAxis();
int nthp = thp_array.size();
/**
@@ -67,35 +69,42 @@ double VFPInjProperties::thp(int table_id,
* 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(table, flo_i, thp_i).value;
const auto flo_i = VFPHelpers<Scalar>::findInterpData(flo, table.getFloAxis());
std::vector<Scalar> bhp_array(nthp);
for (int i = 0; i < nthp; ++i) {
auto thp_i = VFPHelpers<Scalar>::findInterpData(thp_array[i], thp_array);
bhp_array[i] = VFPHelpers<Scalar>::interpolate(table, flo_i, thp_i).value;
}
double retval = detail::findTHP(bhp_array, thp_array, bhp_arg);
return retval;
return VFPHelpers<Scalar>::findTHP(bhp_array, thp_array, bhp_arg);
}
const VFPInjTable& VFPInjProperties::getTable(const int table_id) const {
template<class Scalar>
const VFPInjTable&
VFPInjProperties<Scalar>::getTable(const int table_id) const
{
return detail::getTable(m_tables, table_id);
}
bool VFPInjProperties::hasTable(const int table_id) const {
template<class Scalar>
bool VFPInjProperties<Scalar>::hasTable(const int table_id) const
{
return detail::hasTable(m_tables, table_id);
}
void VFPInjProperties::addTable(const VFPInjTable& new_table) {
template<class Scalar>
void VFPInjProperties<Scalar>::addTable(const VFPInjTable& new_table)
{
this->m_tables.emplace( new_table.getTableNum(), new_table );
}
template<class Scalar>
template <class EvalWell>
EvalWell VFPInjProperties::bhp(const int table_id,
const EvalWell& aqua,
const EvalWell& liquid,
const EvalWell& vapour,
const double& thp) const
EvalWell VFPInjProperties<Scalar>::bhp(const int table_id,
const EvalWell& aqua,
const EvalWell& liquid,
const EvalWell& vapour,
const Scalar thp) const
{
//Get the table
const VFPInjTable& table = detail::getTable(m_tables, table_id);
@@ -106,22 +115,24 @@ EvalWell VFPInjProperties::bhp(const int table_id,
//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(), table.getFloAxis());
auto thp_i = detail::findInterpData( thp, table.getTHPAxis()); // assume constant
const auto flo_i = VFPHelpers<Scalar>::findInterpData(flo.value(), table.getFloAxis());
const auto thp_i = VFPHelpers<Scalar>::findInterpData(thp, table.getTHPAxis()); // assume constant
detail::VFPEvaluation bhp_val = detail::interpolate(table, flo_i, thp_i);
detail::VFPEvaluation bhp_val = VFPHelpers<Scalar>::interpolate(table, flo_i, thp_i);
bhp = bhp_val.dflo * flo;
bhp.setValue(bhp_val.value); // thp is assumed constant i.e.
return bhp;
}
template class VFPInjProperties<double>;
#define INSTANCE(...) \
template __VA_ARGS__ VFPInjProperties::bhp<__VA_ARGS__>(const int, \
const __VA_ARGS__&, \
const __VA_ARGS__&, \
const __VA_ARGS__&, \
const double&) const;
template __VA_ARGS__ VFPInjProperties<double>::bhp<__VA_ARGS__>(const int, \
const __VA_ARGS__&, \
const __VA_ARGS__&, \
const __VA_ARGS__&, \
const double) const;
INSTANCE(DenseAd::Evaluation<double, -1, 4u>)
INSTANCE(DenseAd::Evaluation<double, -1, 5u>)
@@ -139,4 +150,5 @@ INSTANCE(DenseAd::Evaluation<double, 7, 0u>)
INSTANCE(DenseAd::Evaluation<double, 8, 0u>)
INSTANCE(DenseAd::Evaluation<double, 9, 0u>)
INSTANCE(DenseAd::Evaluation<double, 10, 0u>)
} //Namespace Opm

31
opm/simulators/wells/VFPInjProperties.hpp
View File

@@ -30,9 +30,9 @@ namespace Opm {
class VFPInjTable;
template<class Scalar>
class VFPInjProperties {
public:
VFPInjProperties() = default;
/**
* Takes *no* ownership of data.
*/
@@ -55,11 +55,11 @@ public:
* input ADB objects.
*/
template <class EvalWell>
EvalWell bhp(const int table_id,
EvalWell bhp(const int table_id,
const EvalWell& aqua,
const EvalWell& liquid,
const EvalWell& vapour,
const double& thp) const;
const Scalar thp) const;
/**
* Returns the table associated with the ID, or throws an exception if
@@ -75,7 +75,8 @@ public:
/**
* Returns true if no vfp tables are in the current map
*/
bool empty() const {
bool empty() const
{
return m_tables.empty();
}
@@ -90,11 +91,11 @@ public:
* @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;
Scalar bhp(const int table_id,
const Scalar aqua,
const Scalar liquid,
const Scalar vapour,
const Scalar thp) const;
/**
* Linear interpolation of thp as a function of the input parameters
@@ -107,19 +108,17 @@ public:
* @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;
Scalar thp(const int table_id,
const Scalar aqua,
const Scalar liquid,
const Scalar vapour,
const Scalar bhp) const;
protected:
// Map which connects the table number with the table itself
std::map<int, std::reference_wrapper<const VFPInjTable>> m_tables;
};
} //namespace

178
opm/simulators/wells/VFPProdProperties.cpp
View File

@@ -31,21 +31,21 @@
namespace Opm {
double VFPProdProperties::thp(int table_id,
const double& aqua,
const double& liquid,
const double& vapour,
const double& bhp_arg,
const double& alq) const {
template<class Scalar>
Scalar VFPProdProperties<Scalar>::
thp(const int table_id,
const Scalar aqua,
const Scalar liquid,
const Scalar vapour,
const Scalar bhp_arg,
const Scalar alq) const
{
const VFPProdTable& table = detail::getTable(m_tables, table_id);
// Find interpolation variables.
double flo = 0.0;
double wfr = 0.0;
double gfr = 0.0;
Scalar flo = 0.0;
Scalar wfr = 0.0;
Scalar gfr = 0.0;
if (aqua == 0.0 && liquid == 0.0 && vapour == 0.0) {
// All zero, likely at initial state.
// Set FLO variable to minimum to avoid extrapolation.
@@ -67,68 +67,75 @@ double VFPProdProperties::thp(int table_id,
* 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());
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(table, flo_i, thp_i, wfr_i, gfr_i, alq_i).value;
auto flo_i = VFPHelpers<Scalar>::findInterpData( flo, table.getFloAxis());
auto wfr_i = VFPHelpers<Scalar>::findInterpData( wfr, table.getWFRAxis());
auto gfr_i = VFPHelpers<Scalar>::findInterpData( gfr, table.getGFRAxis());
auto alq_i = VFPHelpers<Scalar>::findInterpData( alq, table.getALQAxis());
std::vector<Scalar> bhp_array(nthp);
for (int i = 0; i < nthp; ++i) {
auto thp_i = VFPHelpers<Scalar>::findInterpData(thp_array[i], thp_array);
bhp_array[i] = VFPHelpers<Scalar>::interpolate(table, flo_i, thp_i, wfr_i, gfr_i, alq_i).value;
}
double retval = detail::findTHP(bhp_array, thp_array, bhp_arg);
return retval;
return VFPHelpers<Scalar>::findTHP(bhp_array, thp_array, bhp_arg);
}
double VFPProdProperties::bhp(int table_id,
const double& aqua,
const double& liquid,
const double& vapour,
const double& thp_arg,
const double& alq,
const double& explicit_wfr,
const double& explicit_gfr,
const bool use_expvfp) const {
template<class Scalar>
Scalar VFPProdProperties<Scalar>::
bhp(const int table_id,
const Scalar aqua,
const Scalar liquid,
const Scalar vapour,
const Scalar thp_arg,
const Scalar alq,
const Scalar explicit_wfr,
const Scalar explicit_gfr,
const bool use_expvfp) const
{
const VFPProdTable& table = detail::getTable(m_tables, table_id);
detail::VFPEvaluation retval = detail::bhp(table, aqua, liquid, vapour, thp_arg, alq, explicit_wfr,explicit_gfr, use_expvfp);
detail::VFPEvaluation retval = VFPHelpers<Scalar>::bhp(table, aqua, liquid, vapour,
thp_arg, alq, explicit_wfr,
explicit_gfr, use_expvfp);
return retval.value;
}
const VFPProdTable& VFPProdProperties::getTable(const int table_id) const {
template<class Scalar>
const VFPProdTable&
VFPProdProperties<Scalar>::getTable(const int table_id) const
{
return detail::getTable(m_tables, table_id);
}
bool VFPProdProperties::hasTable(const int table_id) const {
template<class Scalar>
bool VFPProdProperties<Scalar>::hasTable(const int table_id) const
{
return detail::hasTable(m_tables, table_id);
}
std::vector<double>
VFPProdProperties::
bhpwithflo(const std::vector<double>& flos,
template<class Scalar>
std::vector<Scalar>
VFPProdProperties<Scalar>::
bhpwithflo(const std::vector<Scalar>& flos,
const int table_id,
const double wfr,
const double gfr,
const double thp,
const double alq,
const double dp) const
const Scalar wfr,
const Scalar gfr,
const Scalar thp,
const Scalar alq,
const Scalar dp) const
{
// Get the table
const VFPProdTable& table = detail::getTable(m_tables, table_id);
const auto thp_i = detail::findInterpData( thp, table.getTHPAxis()); // assume constant
const auto wfr_i = detail::findInterpData( wfr, table.getWFRAxis());
const auto gfr_i = detail::findInterpData( gfr, table.getGFRAxis());
const auto alq_i = detail::findInterpData( alq, table.getALQAxis()); //assume constant
const auto thp_i = VFPHelpers<Scalar>::findInterpData( thp, table.getTHPAxis()); // assume constant
const auto wfr_i = VFPHelpers<Scalar>::findInterpData( wfr, table.getWFRAxis());
const auto gfr_i = VFPHelpers<Scalar>::findInterpData( gfr, table.getGFRAxis());
const auto alq_i = VFPHelpers<Scalar>::findInterpData( alq, table.getALQAxis()); //assume constant
std::vector<double> bhps(flos.size(), 0.);
std::vector<Scalar> bhps(flos.size(), 0.);
for (std::size_t i = 0; i < flos.size(); ++i) {
// Value of FLO is negative in OPM for producers, but positive in VFP table
const auto flo_i = detail::findInterpData(-flos[i], table.getFloAxis());
const detail::VFPEvaluation bhp_val = detail::interpolate(table, flo_i, thp_i, wfr_i, gfr_i, alq_i);
const auto flo_i = VFPHelpers<Scalar>::findInterpData(-flos[i], table.getFloAxis());
const detail::VFPEvaluation bhp_val = VFPHelpers<Scalar>::interpolate(table, flo_i, thp_i, wfr_i, gfr_i, alq_i);
// TODO: this kind of breaks the conventions for the functions here by putting dp within the function
bhps[i] = bhp_val.value - dp;
@@ -137,37 +144,39 @@ bhpwithflo(const std::vector<double>& flos,
return bhps;
}
double
VFPProdProperties::
template<class Scalar>
Scalar VFPProdProperties<Scalar>::
minimumBHP(const int table_id,
const double thp,
const double wfr,
const double gfr,
const double alq) const
const Scalar thp,
const Scalar wfr,
const Scalar gfr,
const Scalar alq) const
{
// Get the table
const VFPProdTable& table = detail::getTable(m_tables, table_id);
const auto retval = detail::getMinimumBHPCoordinate(table, thp, wfr, gfr, alq);
const auto retval = VFPHelpers<Scalar>::getMinimumBHPCoordinate(table, thp, wfr, gfr, alq);
// returned pair is (flo, bhp)
return retval.second;
}
void VFPProdProperties::addTable(const VFPProdTable& new_table) {
template<class Scalar>
void VFPProdProperties<Scalar>::addTable(const VFPProdTable& new_table)
{
this->m_tables.emplace( new_table.getTableNum(), new_table );
}
template<class Scalar>
template <class EvalWell>
EvalWell VFPProdProperties::bhp(const int table_id,
const EvalWell& aqua,
const EvalWell& liquid,
const EvalWell& vapour,
const double& thp,
const double& alq,
const double& explicit_wfr,
const double& explicit_gfr,
const bool use_expvfp) const
EvalWell VFPProdProperties<Scalar>::
bhp(const int table_id,
const EvalWell& aqua,
const EvalWell& liquid,
const EvalWell& vapour,
const Scalar thp,
const Scalar alq,
const Scalar explicit_wfr,
const Scalar explicit_gfr,
const bool use_expvfp) const
{
//Get the table
const VFPProdTable& table = detail::getTable(m_tables, table_id);
@@ -184,23 +193,32 @@ EvalWell VFPProdProperties::bhp(const int table_id,
//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(), table.getFloAxis());
auto thp_i = detail::findInterpData( thp, table.getTHPAxis()); // assume constant
auto wfr_i = detail::findInterpData( wfr.value(), table.getWFRAxis());
auto gfr_i = detail::findInterpData( gfr.value(), table.getGFRAxis());
auto alq_i = detail::findInterpData( alq, table.getALQAxis()); //assume constant
auto flo_i = VFPHelpers<Scalar>::findInterpData(-flo.value(), table.getFloAxis());
auto thp_i = VFPHelpers<Scalar>::findInterpData( thp, table.getTHPAxis()); // assume constant
auto wfr_i = VFPHelpers<Scalar>::findInterpData( wfr.value(), table.getWFRAxis());
auto gfr_i = VFPHelpers<Scalar>::findInterpData( gfr.value(), table.getGFRAxis());
auto alq_i = VFPHelpers<Scalar>::findInterpData( alq, table.getALQAxis()); //assume constant
detail::VFPEvaluation bhp_val = detail::interpolate(table, flo_i, thp_i, wfr_i, gfr_i, alq_i);
detail::VFPEvaluation bhp_val = VFPHelpers<Scalar>::interpolate(table, flo_i, thp_i, wfr_i,
gfr_i, alq_i);
bhp = (bhp_val.dwfr * wfr) + (bhp_val.dgfr * gfr) - (std::max(0.0, bhp_val.dflo) * flo);
bhp.setValue(bhp_val.value);
return bhp;
}
template class VFPProdProperties<double>;
#define INSTANCE(...) \
template __VA_ARGS__ VFPProdProperties::bhp<__VA_ARGS__>(const int, \
const __VA_ARGS__&, const __VA_ARGS__&, const __VA_ARGS__&, \
const double&, const double&, const double&, const double&, const bool) const;
template __VA_ARGS__ VFPProdProperties<double>::bhp<__VA_ARGS__>(const int, \
const __VA_ARGS__&, \
const __VA_ARGS__&, \
const __VA_ARGS__&, \
const double, \
const double, \
const double, \
const double, \
const bool) const;
INSTANCE(DenseAd::Evaluation<double, -1, 4u>)
INSTANCE(DenseAd::Evaluation<double, -1, 5u>)

70
opm/simulators/wells/VFPProdProperties.hpp
View File

@@ -34,9 +34,9 @@ class VFPProdTable;
* 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.
*/
template<class Scalar>
class VFPProdProperties {
public:
VFPProdProperties() = default;
/**
* Takes *no* ownership of data.
*/
@@ -60,15 +60,15 @@ public:
* input ADB objects.
*/
template <class EvalWell>
EvalWell bhp(const int table_id,
EvalWell bhp(const int table_id,
const EvalWell& aqua,
const EvalWell& liquid,
const EvalWell& vapour,
const double& thp,
const double& alq,
const double& explicit_wfr,
const double& explicit_gfr,
const bool use_expvfp) const;
const Scalar thp,
const Scalar alq,
const Scalar explicit_wfr,
const Scalar explicit_gfr,
const bool use_expvfp) const;
/**
* Linear interpolation of bhp as a function of the input parameters
@@ -82,15 +82,15 @@ public:
* @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 double& explicit_wfr,
const double& explicit_gfr,
const bool use_expvfp) const;
Scalar bhp(const int table_id,
const Scalar aqua,
const Scalar liquid,
const Scalar vapour,
const Scalar thp,
const Scalar alq,
const Scalar explicit_wfr,
const Scalar explicit_gfr,
const bool use_expvfp) const;
/**
* Linear interpolation of thp as a function of the input parameters
@@ -104,12 +104,12 @@ public:
* @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;
Scalar thp(const int table_id,
const Scalar aqua,
const Scalar liquid,
const Scalar vapour,
const Scalar bhp,
const Scalar alq) const;
/**
* Returns the table associated with the ID, or throws an exception if
@@ -125,33 +125,31 @@ public:
/**
* Returns true if no vfp tables are in the current map
*/
bool empty() const {
bool empty() const
{
return m_tables.empty();
}
/**
* Returns minimum bhp for given thp, wfr, gfr and alq
*/
double minimumBHP(const int table_id, const double thp,
const double wfr, const double gfr, const double alq) const;
Scalar minimumBHP(const int table_id, const Scalar thp,
const Scalar wfr, const Scalar gfr, const Scalar alq) const;
protected:
// calculate a group bhp values with a group of flo rate values
std::vector<double> bhpwithflo(const std::vector<double>& flos,
std::vector<Scalar> bhpwithflo(const std::vector<Scalar>& flos,
const int table_id,
const double wfr,
const double gfr,
const double thp,
const double alq,
const double dp) const;
const Scalar wfr,
const Scalar gfr,
const Scalar thp,
const Scalar alq,
const Scalar dp) const;
// Map which connects the table number with the table itself
std::map<int, std::reference_wrapper<const VFPProdTable>> m_tables;
};
} //namespace
} // namespace Opm
#endif /* OPM_AUTODIFF_VFPPRODPROPERTIES_HPP_ */

28
opm/simulators/wells/VFPProperties.hpp
View File

@@ -26,7 +26,6 @@
#include <opm/simulators/wells/VFPHelpers.hpp>
#include <cstddef>
#include <map>
namespace Opm {
@@ -37,6 +36,7 @@ class VFPProdTable;
* A thin wrapper class that holds one VFPProdProperties and one
* VFPInjProperties object.
*/
template<class Scalar>
class VFPProperties {
public:
/**
@@ -48,8 +48,8 @@ public:
VFPProperties(const std::vector<std::reference_wrapper<const VFPInjTable>>& inj_tables,
const std::vector<std::reference_wrapper<const VFPProdTable>>& prod_tables,
const WellState<double>& well_state)
:well_state_(well_state)
const WellState<Scalar>& well_state)
: well_state_(well_state)
{
for (const auto& vfpinj : inj_tables)
this->m_inj.addTable( vfpinj );
@@ -61,18 +61,21 @@ public:
/**
* Returns the VFP properties for injection wells
*/
const VFPInjProperties* getInj() const {
const VFPInjProperties<Scalar>* getInj() const
{
return &m_inj;
}
/**
* Returns the VFP properties for production wells
*/
const VFPProdProperties* getProd() const {
const VFPProdProperties<Scalar>* getProd() const
{
return &m_prod;
}
double getExplicitWFR(const int table_id, const std::size_t well_index) const {
Scalar getExplicitWFR(const int table_id, const std::size_t well_index) const
{
const auto& rates = well_state_.well(well_index).prev_surface_rates;
const auto& pu = well_state_.phaseUsage();
const auto& aqua = pu.phase_used[BlackoilPhases::Aqua]? rates[pu.phase_pos[BlackoilPhases::Aqua]]:0.0;
@@ -82,7 +85,8 @@ public:
return detail::getWFR(table, aqua, liquid, vapour);
}
double getExplicitGFR(const int table_id, const std::size_t well_index) const {
Scalar getExplicitGFR(const int table_id, const std::size_t well_index) const
{
const auto& rates = well_state_.well(well_index).prev_surface_rates;
const auto& pu = well_state_.phaseUsage();
const auto& aqua = pu.phase_used[BlackoilPhases::Aqua]? rates[pu.phase_pos[BlackoilPhases::Aqua]]:0.0;
@@ -93,13 +97,11 @@ public:
}
private:
VFPInjProperties m_inj;
VFPProdProperties m_prod;
const WellState<double>& well_state_;
VFPInjProperties<Scalar> m_inj;
VFPProdProperties<Scalar> m_prod;
const WellState<Scalar>& well_state_;
};
} //Namespace
} // namespace Opm
#endif /* OPM_AUTODIFF_VFPPROPERTIES_HPP_ */

9
opm/simulators/wells/WellBhpThpCalculator.cpp
View File

@@ -919,7 +919,9 @@ isStableSolution(const WellState<Scalar>& well_state,
const auto& table = well_.vfpProperties()->getProd()->getTable(controls.vfp_table_number);
const bool use_vfpexplicit = well_.useVfpExplicit();
detail::VFPEvaluation bhp = detail::bhp(table, aqua, liquid, vapour, thp, well_.getALQ(well_state), wfr, gfr, use_vfpexplicit);
auto bhp = VFPHelpers<double>::bhp(table, aqua, liquid, vapour, thp,
well_.getALQ(well_state), wfr, gfr,
use_vfpexplicit);
if (bhp.dflo >= 0) {
return true;
@@ -964,7 +966,10 @@ estimateStableBhp(const WellState<Scalar>& well_state,
auto bhp_adjusted = [this, &thp, &dp_hydro](const Scalar bhp) {
return bhp - dp_hydro + getVfpBhpAdjustment(bhp, thp);
};
const auto retval = detail::intersectWithIPR(table, thp, wfr, gfr, well_.getALQ(well_state), ipr.first, ipr.second, bhp_adjusted);
const auto retval = VFPHelpers<double>::intersectWithIPR(table, thp, wfr, gfr,
well_.getALQ(well_state),
ipr.first, ipr.second,
bhp_adjusted);
if (retval.has_value()) {
// returned pair is (flo, bhp)
return retval.value().second;

2
opm/simulators/wells/WellGroupHelpers.cpp
View File

@@ -811,7 +811,7 @@ WellGroupHelpers<Scalar>::
computeNetworkPressures(const Network::ExtNetwork& network,
const WellState<Scalar>& well_state,
const GroupState<Scalar>& group_state,
const VFPProdProperties& vfp_prod_props,
const VFPProdProperties<Scalar>& vfp_prod_props,
const Schedule& schedule,
const int report_time_step)
{

4
opm/simulators/wells/WellGroupHelpers.hpp
View File

@@ -37,7 +37,7 @@ template<class Scalar> class GroupState;
namespace Network { class ExtNetwork; }
struct PhaseUsage;
class Schedule;
class VFPProdProperties;
template<class Scalar> class VFPProdProperties;
template<class Scalar> class WellState;
class FieldPropsManager;
@@ -202,7 +202,7 @@ public:
computeNetworkPressures(const Network::ExtNetwork& network,
const WellState<Scalar>& well_state,
const GroupState<Scalar>& group_state,
const VFPProdProperties& vfp_prod_props,
const VFPProdProperties<Scalar>& vfp_prod_props,
const Schedule& schedule,
const int report_time_step);

2
opm/simulators/wells/WellInterfaceGeneric.cpp
View File

@@ -374,7 +374,7 @@ closeCompletions(const WellTestState& wellTestState)
template<class Scalar>
void WellInterfaceGeneric<Scalar>::
setVFPProperties(const VFPProperties* vfp_properties_arg)
setVFPProperties(const VFPProperties<Scalar>* vfp_properties_arg)
{
vfp_properties_ = vfp_properties_arg;
}

8
opm/simulators/wells/WellInterfaceGeneric.hpp
View File

@@ -40,7 +40,7 @@ class ParallelWellInfo;
struct PerforationData;
struct PhaseUsage;
class SummaryState;
class VFPProperties;
template<class Scalar> class VFPProperties;
class WellTestState;
template<class Scalar> class WellState;
template<class Scalar> class SingleWellState;
@@ -94,7 +94,7 @@ public:
void initCompletions();
void closeCompletions(const WellTestState& wellTestState);
void setVFPProperties(const VFPProperties* vfp_properties_arg);
void setVFPProperties(const VFPProperties<Scalar>* vfp_properties_arg);
void setPrevSurfaceRates(WellState<Scalar>& well_state,
const WellState<Scalar>& prev_well_state) const;
void setGuideRate(const GuideRate* guide_rate_arg);
@@ -129,7 +129,7 @@ public:
Scalar gravity() const { return gravity_; }
const VFPProperties* vfpProperties() const { return vfp_properties_; }
const VFPProperties<Scalar>* vfpProperties() const { return vfp_properties_; }
const ParallelWellInfo& parallelWellInfo() const { return parallel_well_info_; }
@@ -358,7 +358,7 @@ protected:
std::vector<Scalar> inj_fc_multiplier_;
Scalar well_efficiency_factor_;
const VFPProperties* vfp_properties_;
const VFPProperties<Scalar>* vfp_properties_;
const GuideRate* guide_rate_;
std::vector<std::string> well_control_log_;

22
tests/test_vfpproperties.cpp
View File

@@ -68,12 +68,12 @@ BOOST_AUTO_TEST_CASE(findInterpData)
double first = 1;
double last = 15;
Opm::detail::InterpData eval0 = Opm::detail::findInterpData(exact, values);
Opm::detail::InterpData eval1 = Opm::detail::findInterpData(interpolate, values);
Opm::detail::InterpData eval2 = Opm::detail::findInterpData(extrapolate_left, values);
Opm::detail::InterpData eval3 = Opm::detail::findInterpData(extrapolate_right, values);
Opm::detail::InterpData eval4 = Opm::detail::findInterpData(first, values);
Opm::detail::InterpData eval5 = Opm::detail::findInterpData(last, values);
auto eval0 = Opm::VFPHelpers<double>::findInterpData(exact, values);
auto eval1 = Opm::VFPHelpers<double>::findInterpData(interpolate, values);
auto eval2 = Opm::VFPHelpers<double>::findInterpData(extrapolate_left, values);
auto eval3 = Opm::VFPHelpers<double>::findInterpData(extrapolate_right, values);
auto eval4 = Opm::VFPHelpers<double>::findInterpData(first, values);
auto eval5 = Opm::VFPHelpers<double>::findInterpData(last, values);
BOOST_CHECK_EQUAL(eval0.ind_[0], 2);
BOOST_CHECK_EQUAL(eval0.ind_[1], 3);
@@ -109,7 +109,7 @@ BOOST_AUTO_TEST_SUITE_END() // HelperTests
* values data is given at
*/
struct TrivialFixture {
typedef Opm::detail::VFPEvaluation VFPEvaluation;
typedef Opm::detail::VFPEvaluation<double> VFPEvaluation;
TrivialFixture() : table_ids(1, 1),
thp_axis{0.0, 1.0},
@@ -208,7 +208,7 @@ struct TrivialFixture {
alq_axis,
data));
properties = std::make_shared<Opm::VFPProdProperties>();
properties = std::make_shared<Opm::VFPProdProperties<double>>();
properties->addTable( *table );
}
@@ -219,7 +219,7 @@ struct TrivialFixture {
std::shared_ptr<Opm::VFPProdProperties> properties;
std::shared_ptr<Opm::VFPProdProperties<double>> properties;
std::shared_ptr<Opm::VFPProdTable> table;
std::vector<int> table_ids;
@@ -612,7 +612,7 @@ VFPPROD \n\
auto deck = parser.parseString(table_str);
bool gaslift_active = false;
Opm::VFPProdTable table(deck["VFPPROD"].front(), gaslift_active, units);
Opm::VFPProdProperties properties;
Opm::VFPProdProperties<double> properties;
properties.addTable( table );
const int n = 5; //Number of points to check per axis
@@ -673,7 +673,7 @@ BOOST_AUTO_TEST_CASE(ParseInterpolateRealisticVFPPROD)
bool gaslift_active = false;
Opm::VFPProdTable table(deck["VFPPROD"].front(), gaslift_active, units);
Opm::VFPProdProperties properties;
Opm::VFPProdProperties<double> properties;
properties.addTable(table);
//Do some rudimentary testing