From 3d86cc366829dd17302c36b132bb877e9b6aaee7 Mon Sep 17 00:00:00 2001 From: Tor Harald Sandve Date: Fri, 9 Sep 2016 11:33:34 +0200 Subject: [PATCH] Add THP support in the denseAD well model Tested on a set of modified SPE1DECK with VFP for injectors and producers --- opm/autodiff/StandardWellsDense.hpp | 229 +++++++++++------- opm/autodiff/VFPHelpers.hpp | 19 +- opm/autodiff/VFPInjProperties.cpp | 29 +++ opm/autodiff/VFPInjProperties.hpp | 9 + opm/autodiff/VFPProdProperties.cpp | 38 ++- opm/autodiff/VFPProdProperties.hpp | 11 + .../WellStateFullyImplicitBlackoil.hpp | 15 +- 7 files changed, 244 insertions(+), 106 deletions(-) diff --git a/opm/autodiff/StandardWellsDense.hpp b/opm/autodiff/StandardWellsDense.hpp index 6e532ca74..934943299 100644 --- a/opm/autodiff/StandardWellsDense.hpp +++ b/opm/autodiff/StandardWellsDense.hpp @@ -260,6 +260,7 @@ namespace Opm { duneB_.mmv(invDCx,Ax); } + // xw = inv(D)*(rw - C*x) void recoverVariable(const BVector& x, BVector& xw) const { BVector resWell = resWell_; duneC_.mmv(x, resWell); @@ -753,33 +754,52 @@ namespace Opm { FluidSystem::waterPvt().inverseFormationVolumeFactor(fs.pvtRegionIndex(), temperature, p_avg); } - //const PhasePresence& perf_cond = (*phase_condition_)[wells().well_cells[perf]]; - if (pu.phase_used[BlackoilPhases::Vapour]) { int gaspos = pu.phase_pos[BlackoilPhases::Vapour] + perf * pu.num_phases; - //if (perf_cond.hasFreeOil()) { - b_perf[gaspos] = FluidSystem::gasPvt().saturatedInverseFormationVolumeFactor(fs.pvtRegionIndex(), temperature, p_avg); - //} -// else { -// double rv = fs.Rv().value; -// b_perf[gaspos] = FluidSystem::gasPvt().inverseFormationVolumeFactor(fs.pvtRegionIndex(), temperature, p_avg, rv); -// } + + if (pu.phase_used[BlackoilPhases::Liquid]) { + int oilpos = pu.phase_pos[BlackoilPhases::Liquid] + perf * pu.num_phases; + const double oilrate = std::abs(xw.wellRates()[oilpos]); //in order to handle negative rates in producers + rvmax_perf[perf] = FluidSystem::gasPvt().saturatedOilVaporizationFactor(fs.pvtRegionIndex(), temperature, p_avg); + if (oilrate > 0) { + const double gasrate = std::abs(xw.wellRates()[gaspos]); + double rv = 0.0; + if (gasrate > 0) { + rv = oilrate / gasrate; + } + rv = std::min(rv, rvmax_perf[perf]); + + b_perf[gaspos] = FluidSystem::gasPvt().inverseFormationVolumeFactor(fs.pvtRegionIndex(), temperature, p_avg, rv); + } + else { + b_perf[gaspos] = FluidSystem::gasPvt().saturatedInverseFormationVolumeFactor(fs.pvtRegionIndex(), temperature, p_avg); + } + + } else { + b_perf[gaspos] = FluidSystem::gasPvt().saturatedInverseFormationVolumeFactor(fs.pvtRegionIndex(), temperature, p_avg); + } } if (pu.phase_used[BlackoilPhases::Liquid]) { int oilpos = pu.phase_pos[BlackoilPhases::Liquid] + perf * pu.num_phases; - //if (perf_cond.hasFreeGas()) { - b_perf[oilpos] = FluidSystem::oilPvt().saturatedInverseFormationVolumeFactor(fs.pvtRegionIndex(), temperature, p_avg); - //} - //else { - // double rs = fs.Rs().value; - // b_perf[oilpos] = FluidSystem::oilPvt().inverseFormationVolumeFactor(fs.pvtRegionIndex(), temperature, p_avg, rs); - //} - } - - if (pu.phase_used[BlackoilPhases::Liquid] && pu.phase_used[BlackoilPhases::Vapour]) { - rsmax_perf[perf] = FluidSystem::oilPvt().saturatedGasDissolutionFactor(fs.pvtRegionIndex(), temperature, p_avg); - rvmax_perf[perf] = FluidSystem::gasPvt().saturatedOilVaporizationFactor(fs.pvtRegionIndex(), temperature, p_avg); + if (pu.phase_used[BlackoilPhases::Vapour]) { + rsmax_perf[perf] = FluidSystem::oilPvt().saturatedGasDissolutionFactor(fs.pvtRegionIndex(), temperature, p_avg); + int gaspos = pu.phase_pos[BlackoilPhases::Vapour] + perf * pu.num_phases; + const double gasrate = std::abs(xw.wellRates()[gaspos]); + if (gasrate > 0) { + const double oilrate = std::abs(xw.wellRates()[oilpos]); + double rs = 0.0; + if (oilrate > 0) { + rs = gasrate / oilrate; + } + rs = std::min(rs, rsmax_perf[perf]); + b_perf[oilpos] = FluidSystem::oilPvt().inverseFormationVolumeFactor(fs.pvtRegionIndex(), temperature, p_avg, rs); + } else { + b_perf[oilpos] = FluidSystem::oilPvt().saturatedInverseFormationVolumeFactor(fs.pvtRegionIndex(), temperature, p_avg); + } + } else { + b_perf[oilpos] = FluidSystem::oilPvt().saturatedInverseFormationVolumeFactor(fs.pvtRegionIndex(), temperature, p_avg); + } } // Surface density. @@ -886,6 +906,7 @@ namespace Opm { // well_state.wellSolutions()[2 * nw + w] = F[Gas]; switch (well_controls_iget_type(wc, current)) { + case THP: case BHP: { //const int sign1 = dxvar_well[w] > 0 ? 1: -1; @@ -908,6 +929,55 @@ namespace Opm { } break; } + + if (well_controls_iget_type(wc, current) == THP) { + + // Do the THP stuff + double aqua = 0.0; + double liquid = 0.0; + double vapour = 0.0; + + const Opm::PhaseUsage& pu = fluid_->phaseUsage(); + + if ((*active_)[ Water ]) { + aqua = well_state.wellRates()[w*np + pu.phase_pos[ Water ] ]; + } + if ((*active_)[ Oil ]) { + liquid = well_state.wellRates()[w*np + pu.phase_pos[ Oil ] ]; + } + if ((*active_)[ Gas ]) { + vapour = well_state.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]; + // pick the density in the top layer + const int perf = wells().well_connpos[w]; + const double rho = well_perforation_densities_[perf]; + + if (well_type == INJECTOR) { + double dp = wellhelpers::computeHydrostaticCorrection( + wells(), w, vfp_properties_->getInj()->getTable(vfp)->getDatumDepth(), + rho, gravity_); + + well_state.bhp()[w] = vfp_properties_->getInj()->bhp(vfp, aqua, liquid, vapour, thp) - dp; + } + else if (well_type == PRODUCER) { + double dp = wellhelpers::computeHydrostaticCorrection( + wells(), w, vfp_properties_->getProd()->getTable(vfp)->getDatumDepth(), + rho, gravity_); + + well_state.bhp()[w] = vfp_properties_->getProd()->bhp(vfp, aqua, liquid, vapour, thp, alq) - dp; + } + else { + OPM_THROW(std::logic_error, "Expected INJECTOR or PRODUCER well"); + } + } + } break; case SURFACE_RATE: @@ -954,65 +1024,7 @@ namespace Opm { break; } - } - - - - - - const Opm::PhaseUsage& pu = fluid_->phaseUsage(); - //Loop over all wells - #pragma omp parallel for schedule(static) - 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 = well_state.wellRates()[w*np + pu.phase_pos[ Water ] ]; - } - if ((*active_)[ Oil ]) { - liquid = well_state.wellRates()[w*np + pu.phase_pos[ Oil ] ]; - } - if ((*active_)[ Gas ]) { - vapour = well_state.wellRates()[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]; - const int perf = wells().well_connpos[w]; - if (well_type == INJECTOR) { - double dp = wellhelpers::computeHydrostaticCorrection( - wells(), w, vfp_properties_->getInj()->getTable(table_id)->getDatumDepth(), - wellPerforationDensities()[perf], gravity_); - - well_state.thp()[w] = vfp_properties_->getInj()->thp(table_id, aqua, liquid, vapour, well_state.bhp()[w] + dp); - } - else if (well_type == PRODUCER) { - double dp = wellhelpers::computeHydrostaticCorrection( - wells(), w, vfp_properties_->getProd()->getTable(table_id)->getDatumDepth(), - wellPerforationDensities()[perf], gravity_); - - well_state.thp()[w] = vfp_properties_->getProd()->thp(table_id, aqua, liquid, vapour, well_state.bhp()[w] + dp, alq); - } - else { - OPM_THROW(std::logic_error, "Expected INJECTOR or PRODUCER well"); - } - - //Assume only one THP control specified for each well - break; - } - } - } + } } } @@ -1102,19 +1114,22 @@ namespace Opm { //Set *BHP* target by calculating bhp from THP const WellType& well_type = wells().type[w]; + + // pick the density in the top layer const int perf = wells().well_connpos[w]; + const double rho = well_perforation_densities_[perf]; if (well_type == INJECTOR) { double dp = wellhelpers::computeHydrostaticCorrection( wells(), w, vfp_properties_->getInj()->getTable(vfp)->getDatumDepth(), - wellPerforationDensities()[perf], gravity_); + rho, gravity_); xw.bhp()[w] = vfp_properties_->getInj()->bhp(vfp, aqua, liquid, vapour, thp) - dp; } else if (well_type == PRODUCER) { double dp = wellhelpers::computeHydrostaticCorrection( wells(), w, vfp_properties_->getProd()->getTable(vfp)->getDatumDepth(), - wellPerforationDensities()[perf], gravity_); + rho, gravity_); xw.bhp()[w] = vfp_properties_->getProd()->bhp(vfp, aqua, liquid, vapour, thp, alq) - dp; } @@ -1173,6 +1188,7 @@ namespace Opm { } } switch (well_controls_iget_type(wc, current)) { + case THP: case BHP: { const WellType& well_type = wells().type[w]; @@ -1203,11 +1219,11 @@ namespace Opm { tot_well_rate += g[p] * xw.wellRates()[np*w + p]; } if(std::abs(tot_well_rate) > 0) { - xw.wellSolutions()[nw + w] = g[Water] * xw.wellRates()[np*w + Water] / tot_well_rate; //wells->comp_frac[np*w + Water]; // Water; - xw.wellSolutions()[2*nw + w] = g[Gas] * xw.wellRates()[np*w + Gas] / tot_well_rate ; //wells->comp_frac[np*w + Gas]; //Gas + xw.wellSolutions()[nw + w] = g[Water] * xw.wellRates()[np*w + Water] / tot_well_rate; + xw.wellSolutions()[2*nw + w] = g[Gas] * xw.wellRates()[np*w + Gas] / tot_well_rate ; } else { - //xw.wellSolutions()[nw + w] = wells().comp_frac[np*w + Water]; - //xw.wellSolutions()[2 * nw + w] = wells().comp_frac[np*w + Gas]; + xw.wellSolutions()[nw + w] = wells().comp_frac[np*w + Water]; + xw.wellSolutions()[2 * nw + w] = wells().comp_frac[np*w + Gas]; } } } @@ -1382,6 +1398,43 @@ namespace Opm { const double target_rate = well_controls_get_current_target(wc); bhp.value = target_rate; return bhp; + } else if (well_controls_get_current_type(wc) == THP) { + const int control = well_controls_get_current(wc); + const double thp = well_controls_get_current_target(wc); + const double alq = well_controls_iget_alq(wc, control); + const int table_id = well_controls_iget_vfp(wc, control); + EvalWell aqua = 0.0; + EvalWell liquid = 0.0; + EvalWell vapour = 0.0; + EvalWell bhp = 0.0; + double vfp_ref_depth = 0.0; + + const Opm::PhaseUsage& pu = fluid_->phaseUsage(); + + if ((*active_)[ Water ]) { + aqua = getQs(wellIdx, pu.phase_pos[ Water]); + } + if ((*active_)[ Oil ]) { + liquid = getQs(wellIdx, pu.phase_pos[ Oil ]); + } + if ((*active_)[ Gas ]) { + vapour = getQs(wellIdx, pu.phase_pos[ Gas ]); + } + if (wells().type[wellIdx] == INJECTOR) { + bhp = vfp_properties_->getInj()->bhp(table_id, aqua, liquid, vapour, thp); + vfp_ref_depth = vfp_properties_->getInj()->getTable(table_id)->getDatumDepth(); + } else { + bhp = vfp_properties_->getProd()->bhp(table_id, aqua, liquid, vapour, thp, alq); + vfp_ref_depth = vfp_properties_->getProd()->getTable(table_id)->getDatumDepth(); + } + + // pick the density in the top layer + const int perf = wells().well_connpos[wellIdx]; + const double rho = well_perforation_densities_[perf]; + const double dp = wellhelpers::computeHydrostaticCorrection(wells(), wellIdx, vfp_ref_depth, rho, gravity_); + bhp -= dp; + return bhp; + } return wellVariables_[wellIdx]; } @@ -1397,7 +1450,7 @@ namespace Opm { if (comp_frac == 0.0) return qs; - if (well_controls_get_current_type(wc) == BHP) { + if (well_controls_get_current_type(wc) == BHP || well_controls_get_current_type(wc) == THP) { return wellVariables_[wellIdx]; } qs.value = target_rate; @@ -1405,7 +1458,7 @@ namespace Opm { } // Producers - if (well_controls_get_current_type(wc) == BHP) { + if (well_controls_get_current_type(wc) == BHP || well_controls_get_current_type(wc) == THP ) { return wellVariables_[wellIdx] * wellVolumeFractionScaled(wellIdx,phaseIdx); } if (well_controls_get_current_type(wc) == SURFACE_RATE) { diff --git a/opm/autodiff/VFPHelpers.hpp b/opm/autodiff/VFPHelpers.hpp index 6220b71de..6760cb7f2 100644 --- a/opm/autodiff/VFPHelpers.hpp +++ b/opm/autodiff/VFPHelpers.hpp @@ -25,7 +25,8 @@ #include #include #include - +#include +#include /** * This file contains a set of helper functions used by VFPProd / VFPInj. @@ -35,14 +36,7 @@ namespace detail { typedef AutoDiffBlock ADB; - - - - - - - - +typedef DenseAd::Evaluation EvalWell; /** @@ -52,7 +46,12 @@ inline double zeroIfNan(const double& value) { return (std::isnan(value)) ? 0.0 : value; } - +/** + * Returns zero if input value is NaN + */ +inline double zeroIfNan(const EvalWell& value) { + return (std::isnan(value.value)) ? 0.0 : value.value; +} diff --git a/opm/autodiff/VFPInjProperties.cpp b/opm/autodiff/VFPInjProperties.cpp index 7027c9957..e2316bc17 100644 --- a/opm/autodiff/VFPInjProperties.cpp +++ b/opm/autodiff/VFPInjProperties.cpp @@ -89,7 +89,36 @@ VFPInjProperties::ADB VFPInjProperties::bhp(const std::vector& table_id, +VFPInjProperties::EvalWell VFPInjProperties::bhp(const int table_id, + const EvalWell& aqua, + const EvalWell& liquid, + const EvalWell& vapour, + const double& thp) const { + //Get the table + const VFPInjTable* table = detail::getTable(m_tables, table_id); + EvalWell bhp = 0.0; + + //Find interpolation variables + EvalWell flo = detail::getFlo(aqua, liquid, vapour, table->getFloType()); + + //Compute the BHP for each well independently + 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, table->getFloAxis()); + auto thp_i = detail::findInterpData( thp, table->getTHPAxis()); // assume constant + + detail::VFPEvaluation bhp_val = detail::interpolate(table->getTable(), flo_i, thp_i); + + bhp = bhp_val.dflo * flo; + bhp.value = bhp_val.value; // thp is assumed constant i.e. + } + else { + bhp.value = -1e100; //Signal that this value has not been calculated properly, due to "missing" table + } + return bhp; +} diff --git a/opm/autodiff/VFPInjProperties.hpp b/opm/autodiff/VFPInjProperties.hpp index 7a6e4c71e..20c79af19 100644 --- a/opm/autodiff/VFPInjProperties.hpp +++ b/opm/autodiff/VFPInjProperties.hpp @@ -23,6 +23,8 @@ #include #include #include +#include +#include #include #include @@ -91,6 +93,13 @@ public: const ADB& vapour, const ADB& thp) const; + typedef DenseAd::Evaluation EvalWell; + EvalWell bhp(const int table_id, + const EvalWell& aqua, + const EvalWell& liquid, + const EvalWell& vapour, + const double& thp) const; + /** * Linear interpolation of bhp as a function of the input parameters * @param table_id Table number to use diff --git a/opm/autodiff/VFPProdProperties.cpp b/opm/autodiff/VFPProdProperties.cpp index 8da911c72..c2bec752f 100644 --- a/opm/autodiff/VFPProdProperties.cpp +++ b/opm/autodiff/VFPProdProperties.cpp @@ -24,7 +24,8 @@ #include #include #include - +#include +#include #include @@ -74,7 +75,42 @@ VFPProdProperties::ADB VFPProdProperties::bhp(const std::vector& table_id, } +VFPProdProperties::EvalWell VFPProdProperties::bhp(const int table_id, + const EvalWell& aqua, + const EvalWell& liquid, + const EvalWell& vapour, + const double& thp, + const double& alq) const { + //Get the table + const VFPProdTable* table = detail::getTable(m_tables, table_id); + EvalWell bhp = 0.0; + + //Find interpolation variables + EvalWell flo = detail::getFlo(aqua, liquid, vapour, table->getFloType()); + EvalWell wfr = detail::getWFR(aqua, liquid, vapour, table->getWFRType()); + EvalWell gfr = detail::getGFR(aqua, liquid, vapour, table->getGFRType()); + + //Compute the BHP for each well independently + 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, 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 + + detail::VFPEvaluation bhp_val = detail::interpolate(table->getTable(), flo_i, thp_i, wfr_i, gfr_i, alq_i); + + bhp = (bhp_val.dwfr * wfr) + (bhp_val.dgfr * gfr) - (bhp_val.dflo * flo); + bhp.value = bhp_val.value; + } + else { + bhp.value = -1e100; //Signal that this value has not been calculated properly, due to "missing" table + } + return bhp; +} diff --git a/opm/autodiff/VFPProdProperties.hpp b/opm/autodiff/VFPProdProperties.hpp index 76580ecdf..33a5bbba7 100644 --- a/opm/autodiff/VFPProdProperties.hpp +++ b/opm/autodiff/VFPProdProperties.hpp @@ -23,6 +23,8 @@ #include #include #include +#include +#include #include #include @@ -99,6 +101,15 @@ public: const ADB& thp, const ADB& alq) const; + + typedef DenseAd::Evaluation EvalWell; + EvalWell bhp(const int table_id, + const EvalWell& aqua, + const EvalWell& liquid, + const EvalWell& vapour, + const double& thp, + const double& alq) const; + /** * Linear interpolation of bhp as a function of the input parameters * @param table_id Table number to use diff --git a/opm/autodiff/WellStateFullyImplicitBlackoil.hpp b/opm/autodiff/WellStateFullyImplicitBlackoil.hpp index 6de120ba8..859f22ce4 100644 --- a/opm/autodiff/WellStateFullyImplicitBlackoil.hpp +++ b/opm/autodiff/WellStateFullyImplicitBlackoil.hpp @@ -119,6 +119,7 @@ namespace Opm const WellType& well_type = wells->type[w]; switch (well_controls_iget_type(wc, current)) { + case THP: // Intentional fall-through case BHP: { if (well_type == INJECTOR) { @@ -148,13 +149,13 @@ namespace Opm total_rates += g[p] * wellRates()[np*w + p]; } - //if(std::abs(total_rates) > 0) { - // wellSolutions()[nw + w] = g[Water] * wellRates()[np*w + Water] / total_rates; //wells->comp_frac[np*w + Water]; // Water; - // wellSolutions()[2*nw + w] = g[Gas] * wellRates()[np*w + Gas] / total_rates ; //wells->comp_frac[np*w + Gas]; //Gas - //} else { - wellSolutions()[nw + w] = wells->comp_frac[np*w + Water]; - wellSolutions()[2*nw + w] = wells->comp_frac[np*w + Gas]; - //} + if(std::abs(total_rates) > 0) { + wellSolutions()[nw + w] = g[Water] * wellRates()[np*w + Water] / total_rates; //wells->comp_frac[np*w + Water]; // Water; + wellSolutions()[2*nw + w] = g[Gas] * wellRates()[np*w + Gas] / total_rates ; //wells->comp_frac[np*w + Gas]; //Gas + } else { + wellSolutions()[nw + w] = wells->comp_frac[np*w + Water]; + wellSolutions()[2*nw + w] = wells->comp_frac[np*w + Gas]; + } }