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

Fixed sign issue with vfpprod->bhp(flo)

Browse Source
This commit is contained in:
babrodtk 2015-08-18 14:53:36 +02:00
parent 00b63f303f
commit b75baac8fc
4 changed files with 47 additions and 20 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
opm/autodiff/BlackoilModelBase.hpp
View File

@ -374,7 +374,7 @@ namespace Opm {
bool getWellConvergence(const int iteration);
bool isVFPActive(const WellState& well_state) const;
bool isVFPActive() const;
std::vector<ADB>
computePressures(const ADB& po,

24
opm/autodiff/BlackoilModelBase_impl.hpp
View File

@ -777,7 +777,7 @@ namespace detail {
// 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(well_state)) {
if (isVFPActive()) {
SolutionState state = asImpl().variableState(reservoir_state, well_state);
SolutionState state0 = state;
asImpl().makeConstantState(state0);
@ -927,7 +927,6 @@ namespace detail {
template <class Grid, class Implementation>
void
BlackoilModelBase<Grid, Implementation>::computeWellFlux(const SolutionState& state,
@ -1191,25 +1190,19 @@ 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.
if (well_perforation_densities.size() > w) {
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;
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;
}
else {
return 0.0;
}
return dp;
}
} //Namespace
template <class Grid, class Implementation>
bool BlackoilModelBase<Grid, Implementation>::isVFPActive(const WellState& xw) const
bool BlackoilModelBase<Grid, Implementation>::isVFPActive() const
{
if( ! wellsActive() ) {
return false;
@ -1576,6 +1569,7 @@ namespace detail {
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;

7
opm/autodiff/VFPHelpers.hpp
View File

@ -208,6 +208,8 @@ struct InterpData {
inline InterpData findInterpData(const double& value, const std::vector<double>& values) {
InterpData retval;
const double abs_value = std::abs(value);
//If we only have one value in our vector, return that
if (values.size() == 1) {
retval.ind_[0] = 0;
@ -220,7 +222,7 @@ inline InterpData findInterpData(const double& value, const std::vector<double>&
//First element greater than or equal to value
//Start with the second element, so that floor_iter does not go out of range
//Don't access out-of-range, therefore values.end()-1
auto ceil_iter = std::lower_bound(values.begin()+1, values.end()-1, value);
auto ceil_iter = std::lower_bound(values.begin()+1, values.end()-1, abs_value);
//Find last element smaller than range
auto floor_iter = ceil_iter-1;
@ -235,7 +237,7 @@ inline InterpData findInterpData(const double& value, const std::vector<double>&
//Possible source for floating point error here if value and floor are large,
//but very close to each other
retval.inv_dist_ = 1.0 / dist;
retval.factor_ = (value-*floor_iter) * retval.inv_dist_;
retval.factor_ = (abs_value-*floor_iter) * retval.inv_dist_;
}
else {
retval.inv_dist_ = 0.0;
@ -519,7 +521,6 @@ inline adb_like bhp(const VFPProdTable* table,
auto gfr_i = detail::findInterpData(gfr, table->getGFRAxis());
auto alq_i = detail::findInterpData(alq, table->getALQAxis());
//Then perform the interpolation itself
detail::adb_like retval = detail::interpolate(table->getTable(), flo_i, thp_i, wfr_i, gfr_i, alq_i);
return retval;

34
opm/autodiff/VFPProdProperties.cpp
View File

@ -130,6 +130,38 @@ VFPProdProperties::ADB VFPProdProperties::bhp(const std::vector<int>& table_id,
detail::adb_like bhp_val = detail::interpolate(table->getTable(), flo_i, thp_i, wfr_i, gfr_i, alq_i);
/*
static const int N=40;
std::cout << "bhp=" << bhp_val.value << ";" << std::endl;
std::cout << "flo=" << flo.value()[i] << ";" << std::endl;
std::cout << "thp=" << thp.value()[i] << ";" << std::endl;
std::cout << "wfr=" << wfr.value()[i] << ";" << std::endl;
std::cout << "gfr=" << gfr.value()[i] << ";" << std::endl;
std::cout << "alq=" << alq.value()[i] << ";" << std::endl;
std::cout << "bhp_vfp=[" << std::endl;
for (int j=0; j<N; ++j) {
const double start = table->getFloAxis().front();
const double end = table->getFloAxis().back();
const double dist = end - start;
double flo_d = start + (j/static_cast<double>(N-1)) * dist;
auto flo_i = detail::findInterpData(flo_d, table->getFloAxis());
detail::adb_like bhp_val = detail::interpolate(table->getTable(), flo_i, thp_i, wfr_i, gfr_i, alq_i);
std::cout << bhp_val.value << ",";
}
std::cout << "];" << std::endl;
std::cout << "flo_vfp=[" << std::endl;
for (int j=0; j<N; ++j) {
const double start = table->getFloAxis().front();
const double end = table->getFloAxis().back();
const double dist = end - start;
double flo_d = start + (j/static_cast<double>(N-1)) * dist;
std::cout << flo_d << ",";
}
std::cout << "];" << std::endl;
*/
value[i] = bhp_val.value;
dthp[i] = bhp_val.dthp;
dwfr[i] = bhp_val.dwfr;
@ -170,7 +202,7 @@ VFPProdProperties::ADB VFPProdProperties::bhp(const std::vector<int>& table_id,
jacs[block] += dalq_diag * alq.derivative()[block];
}
if (!flo.derivative().empty()) {
jacs[block] += dflo_diag * flo.derivative()[block];
jacs[block] -= dflo_diag * flo.derivative()[block];
}
}