diff --git a/opm/autodiff/AutoDiffHelpers.hpp b/opm/autodiff/AutoDiffHelpers.hpp index 2c7d5b54a..ce486f20d 100644 --- a/opm/autodiff/AutoDiffHelpers.hpp +++ b/opm/autodiff/AutoDiffHelpers.hpp @@ -365,7 +365,7 @@ spdiag(const AutoDiffBlock::V& d) public: typedef AutoDiffBlock ADB; - enum CriterionForLeftElement { GreaterEqualZero, GreaterZero, Zero, NotEqualZero, LessZero, LessEqualZero, NotNaN }; + enum CriterionForLeftElement { GreaterEqualZero, GreaterZero, Zero, NotEqualZero, LessZero, LessEqualZero, NotNaN, NotNaNInf }; Selector(const typename ADB::V& selection_basis, CriterionForLeftElement crit = GreaterEqualZero) @@ -400,6 +400,9 @@ spdiag(const AutoDiffBlock::V& d) case NotNaN: chooseleft = !isnan(selection_basis[i]); break; + case NotNaNInf: + chooseleft = !isnan(selection_basis[i]) && !std::isinf(selection_basis[i]); + break; default: OPM_THROW(std::logic_error, "No such criterion: " << crit); } diff --git a/opm/autodiff/BlackoilWellModel_impl.hpp b/opm/autodiff/BlackoilWellModel_impl.hpp index 3ea869fcc..34b30fcd9 100644 --- a/opm/autodiff/BlackoilWellModel_impl.hpp +++ b/opm/autodiff/BlackoilWellModel_impl.hpp @@ -844,8 +844,6 @@ namespace Opm { // we simply return. if( !wellsActive() ) return ; -#if HAVE_OPENMP -#endif // HAVE_OPENMP wellhelpers::WellSwitchingLogger logger; for (const auto& well : well_container_) { diff --git a/opm/autodiff/StandardWell.hpp b/opm/autodiff/StandardWell.hpp index 234dd6a4e..0dad16f20 100644 --- a/opm/autodiff/StandardWell.hpp +++ b/opm/autodiff/StandardWell.hpp @@ -331,7 +331,7 @@ namespace Opm template ValueType calculateBhpFromThp(const std::vector& rates, const int control_index) const; - double calculateThpFromBhp(const std::vector& rates, const int control_index, const double bhp) const; + double calculateThpFromBhp(const std::vector& rates, const double bhp) const; // get the mobility for specific perforation void getMobility(const Simulator& ebosSimulator, diff --git a/opm/autodiff/StandardWell_impl.hpp b/opm/autodiff/StandardWell_impl.hpp index 384c3ebaa..6ec57c871 100644 --- a/opm/autodiff/StandardWell_impl.hpp +++ b/opm/autodiff/StandardWell_impl.hpp @@ -1084,43 +1084,35 @@ namespace Opm StandardWell:: updateThp(WellState& well_state) const { - // for the wells having a THP constaint, we should update their thp value - // If it is under THP control, it will be set to be the target value. - // TODO: a better standard is probably whether we have the table to calculate the THP value - // TODO: it is something we need to check the output to decide. - const WellControls* wc = well_controls_; - // TODO: we should only maintain one current control either from the well_state or from well_controls struct. - // Either one can be more favored depending on the final strategy for the initilzation of the well control - const int nwc = well_controls_get_num(wc); - // Looping over all controls until we find a THP constraint - for (int ctrl_index = 0; ctrl_index < nwc; ++ctrl_index) { - if (well_controls_iget_type(wc, ctrl_index) == THP) { - // the current control - const int current = well_state.currentControls()[index_of_well_]; - // if well under THP control at the moment - if (current == ctrl_index) { - const double thp_target = well_controls_iget_target(wc, current); - well_state.thp()[index_of_well_] = thp_target; - } else { // otherwise we calculate the thp from the bhp value - const Opm::PhaseUsage& pu = phaseUsage(); - std::vector rates(3, 0.0); + // When there is no vaild VFP table provided, we set the thp to be zero. + if (!this->isVFPActive()) { + well_state.thp()[index_of_well_] = 0.; + return; + } - if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) { - rates[ Water ] = well_state.wellRates()[index_of_well_ * number_of_phases_ + pu.phase_pos[ Water ] ]; - } - if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) { - rates[ Oil ] = well_state.wellRates()[index_of_well_ * number_of_phases_ + pu.phase_pos[ Oil ] ]; - } - if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) { - rates[ Gas ] = well_state.wellRates()[index_of_well_ * number_of_phases_ + pu.phase_pos[ Gas ] ]; - } + // avaiable VFP table is provided, we should update the thp value - const double bhp = well_state.bhp()[index_of_well_]; + // if the well is under THP control, we should use its target value + if (well_controls_get_current_type(well_controls_) == THP) { + well_state.thp()[index_of_well_] = well_controls_get_current_target(well_controls_); + } else { + // the well is under other control types, we calculate the thp based on bhp and rates + std::vector rates(3, 0.0); - well_state.thp()[index_of_well_] = calculateThpFromBhp(rates, ctrl_index, bhp); - } - break; + const Opm::PhaseUsage& pu = phaseUsage(); + if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) { + rates[ Water ] = well_state.wellRates()[index_of_well_ * number_of_phases_ + pu.phase_pos[ Water ] ]; } + if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) { + rates[ Oil ] = well_state.wellRates()[index_of_well_ * number_of_phases_ + pu.phase_pos[ Oil ] ]; + } + if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) { + rates[ Gas ] = well_state.wellRates()[index_of_well_ * number_of_phases_ + pu.phase_pos[ Gas ] ]; + } + + const double bhp = well_state.bhp()[index_of_well_]; + + well_state.thp()[index_of_well_] = calculateThpFromBhp(rates, bhp); } } @@ -2096,7 +2088,6 @@ namespace Opm double StandardWell:: calculateThpFromBhp(const std::vector& rates, - const int control_index, const double bhp) const { assert(int(rates.size()) == 3); // the vfp related only supports three phases now. @@ -2105,32 +2096,30 @@ namespace Opm const double liquid = rates[Oil]; const double vapour = rates[Gas]; - const int vfp = well_controls_iget_vfp(well_controls_, control_index); - const double& alq = well_controls_iget_alq(well_controls_, control_index); - // pick the density in the top layer const double rho = perf_densities_[0]; double thp = 0.0; if (well_type_ == INJECTOR) { - const double vfp_ref_depth = vfp_properties_->getInj()->getTable(vfp)->getDatumDepth(); - + const int table_id = well_ecl_->getInjectionProperties(current_step_).VFPTableNumber; + const double vfp_ref_depth = vfp_properties_->getInj()->getTable(table_id)->getDatumDepth(); const double dp = wellhelpers::computeHydrostaticCorrection(ref_depth_, vfp_ref_depth, rho, gravity_); - thp = vfp_properties_->getInj()->thp(vfp, aqua, liquid, vapour, bhp + dp); - } - else if (well_type_ == PRODUCER) { - const double vfp_ref_depth = vfp_properties_->getProd()->getTable(vfp)->getDatumDepth(); + thp = vfp_properties_->getInj()->thp(table_id, aqua, liquid, vapour, bhp + dp); + } + else if (well_type_ == PRODUCER) { + const int table_id = well_ecl_->getProductionProperties(current_step_).VFPTableNumber; + const double alq = well_ecl_->getProductionProperties(current_step_).ALQValue; + const double vfp_ref_depth = vfp_properties_->getProd()->getTable(table_id)->getDatumDepth(); + const double dp = wellhelpers::computeHydrostaticCorrection(ref_depth_, vfp_ref_depth, rho, gravity_); - const double dp = wellhelpers::computeHydrostaticCorrection(ref_depth_, vfp_ref_depth, rho, gravity_); + thp = vfp_properties_->getProd()->thp(table_id, aqua, liquid, vapour, bhp + dp, alq); + } + else { + OPM_THROW(std::logic_error, "Expected INJECTOR or PRODUCER well"); + } - thp = vfp_properties_->getProd()->thp(vfp, aqua, liquid, vapour, bhp + dp, alq); - } - else { - OPM_THROW(std::logic_error, "Expected INJECTOR or PRODUCER well"); - } - - return thp; + return thp; } diff --git a/opm/autodiff/VFPHelpers.hpp b/opm/autodiff/VFPHelpers.hpp index 80f4b1c38..81a4916f1 100644 --- a/opm/autodiff/VFPHelpers.hpp +++ b/opm/autodiff/VFPHelpers.hpp @@ -21,6 +21,7 @@ #ifndef OPM_AUTODIFF_VFPHELPERS_HPP_ #define OPM_AUTODIFF_VFPHELPERS_HPP_ +#include #include #include @@ -37,19 +38,31 @@ namespace detail { /** - * Returns zero if input value is NaN + * Returns zero if input value is NaN of INF */ -inline double zeroIfNan(const double& value) { - return (std::isnan(value)) ? 0.0 : value; +inline double zeroIfNanInf(const double& value) { + const bool nan_or_inf = std::isnan(value) || std::isinf(value); + + if (nan_or_inf) { + OpmLog::warning("NAN_OR_INF_VFP", "NAN or INF value encountered during VFP calculation, the value is set to zero"); + } + + return nan_or_inf ? 0.0 : value; } /** - * Returns zero if input value is NaN + * Returns zero if input value is NaN or INF */ template -inline EvalWell zeroIfNan(const EvalWell& value) { - return (std::isnan(value.value())) ? 0.0 : value; +inline EvalWell zeroIfNanInf(const EvalWell& value) { + const bool nan_or_inf = std::isnan(value.value()) || std::isinf(value.value()); + + if (nan_or_inf) { + OpmLog::warning("NAN_OR_INF_VFP_EVAL", "NAN or INF Evalution encountered during VFP calculation, the Evalution is set to zero"); + } + + return nan_or_inf ? 0.0 : value; } @@ -120,14 +133,14 @@ static T getWFR(const T& aqua, const T& liquid, const T& vapour, case VFPProdTable::WFR_WOR: { //Water-oil ratio = water / oil T wor = aqua / liquid; - return zeroIfNan(wor); + return zeroIfNanInf(wor); } case VFPProdTable::WFR_WCT: //Water cut = water / (water + oil) - return zeroIfNan(aqua / (aqua + liquid)); + return zeroIfNanInf(aqua / (aqua + liquid)); case VFPProdTable::WFR_WGR: //Water-gas ratio = water / gas - return zeroIfNan(aqua / vapour); + return zeroIfNanInf(aqua / vapour); case VFPProdTable::WFR_INVALID: //Intentional fall-through default: OPM_THROW(std::logic_error, "Invalid WFR_TYPE: '" << type << "'"); @@ -149,13 +162,13 @@ static T getGFR(const T& aqua, const T& liquid, const T& vapour, switch(type) { case VFPProdTable::GFR_GOR: // Gas-oil ratio = gas / oil - return zeroIfNan(vapour / liquid); + return zeroIfNanInf(vapour / liquid); case VFPProdTable::GFR_GLR: // Gas-liquid ratio = gas / (oil + water) - return zeroIfNan(vapour / (liquid + aqua)); + return zeroIfNanInf(vapour / (liquid + aqua)); case VFPProdTable::GFR_OGR: // Oil-gas ratio = oil / gas - return zeroIfNan(liquid / vapour); + return zeroIfNanInf(liquid / vapour); case VFPProdTable::GFR_INVALID: //Intentional fall-through default: OPM_THROW(std::logic_error, "Invalid GFR_TYPE: '" << type << "'"); @@ -537,13 +550,22 @@ template const T* getTable(const std::map tables, int table_id) { auto entry = tables.find(table_id); if (entry == tables.end()) { - OPM_THROW(std::invalid_argument, "Nonexistent table " << table_id << " referenced."); + OPM_THROW(std::invalid_argument, "Nonexistent VFP table " << table_id << " referenced."); } else { return entry->second; } } +/** + * Check whether we have a table with the table number + */ +template +bool hasTable(const std::map tables, int table_id) { + const auto entry = tables.find(table_id); + return (entry != tables.end() ); +} + /** * Returns the type variable for FLO/GFR/WFR for production tables diff --git a/opm/autodiff/VFPHelpersLegacy.hpp b/opm/autodiff/VFPHelpersLegacy.hpp index bb8a9f3d7..b9583d1a8 100644 --- a/opm/autodiff/VFPHelpersLegacy.hpp +++ b/opm/autodiff/VFPHelpersLegacy.hpp @@ -37,15 +37,15 @@ typedef AutoDiffBlock ADB; /** - * Returns zero for every entry in the ADB which is NaN - */ -inline ADB zeroIfNan(const ADB& values) { - Selector not_nan_selector(values.value(), Selector::NotNaN); + * * Returns zero for every entry in the ADB which is NaN or INF + * */ +inline ADB zeroIfNanInf(const ADB& values) { + Selector not_nan_inf_selector(values.value(), Selector::NotNaNInf); 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); + ADB retval = not_nan_inf_selector.select(values, zero); return retval; } diff --git a/opm/autodiff/VFPInjProperties.cpp b/opm/autodiff/VFPInjProperties.cpp index fbc669370..8d364e8dd 100644 --- a/opm/autodiff/VFPInjProperties.cpp +++ b/opm/autodiff/VFPInjProperties.cpp @@ -91,9 +91,12 @@ double VFPInjProperties::thp(int table_id, return retval; } - const VFPInjTable* VFPInjProperties::getTable(const int table_id) const { return detail::getTable(m_tables, table_id); } +bool VFPInjProperties::hasTable(const int table_id) const { + return detail::hasTable(m_tables, table_id); +} + } //Namespace Opm diff --git a/opm/autodiff/VFPInjProperties.hpp b/opm/autodiff/VFPInjProperties.hpp index ef319442f..044d7a90b 100644 --- a/opm/autodiff/VFPInjProperties.hpp +++ b/opm/autodiff/VFPInjProperties.hpp @@ -109,6 +109,11 @@ public: */ const VFPInjTable* getTable(const int table_id) const; + /** + * Check whether there is table associated with ID + */ + bool hasTable(const int table_id) const; + /** * Returns true if no vfp tables are in the current map */ diff --git a/opm/autodiff/VFPProdProperties.cpp b/opm/autodiff/VFPProdProperties.cpp index 924bc60cc..5e17f2dee 100644 --- a/opm/autodiff/VFPProdProperties.cpp +++ b/opm/autodiff/VFPProdProperties.cpp @@ -105,5 +105,9 @@ const VFPProdTable* VFPProdProperties::getTable(const int table_id) const { return detail::getTable(m_tables, table_id); } +bool VFPProdProperties::hasTable(const int table_id) const { + return detail::hasTable(m_tables, table_id); +} + } diff --git a/opm/autodiff/VFPProdProperties.hpp b/opm/autodiff/VFPProdProperties.hpp index 966ade863..b68218a84 100644 --- a/opm/autodiff/VFPProdProperties.hpp +++ b/opm/autodiff/VFPProdProperties.hpp @@ -158,6 +158,11 @@ public: */ const VFPProdTable* getTable(const int table_id) const; + /** + * Check whether there is table associated with ID + */ + bool hasTable(const int table_id) const; + /** * Returns true if no vfp tables are in the current map */ diff --git a/opm/autodiff/WellInterface.hpp b/opm/autodiff/WellInterface.hpp index e2b59e4cf..f2ed3975d 100644 --- a/opm/autodiff/WellInterface.hpp +++ b/opm/autodiff/WellInterface.hpp @@ -334,6 +334,9 @@ namespace Opm double scalingFactor(const int comp_idx) const; + // whether a well is specified with a non-zero and valid VFP table number + bool isVFPActive() const; + void wellTestingEconomic(Simulator& simulator, const std::vector& B_avg, const double simulation_time, const int report_step, const bool terminal_output, const WellState& well_state, WellTestState& welltest_state); diff --git a/opm/autodiff/WellInterface_impl.hpp b/opm/autodiff/WellInterface_impl.hpp index 22588c3fa..08e5d46d8 100644 --- a/opm/autodiff/WellInterface_impl.hpp +++ b/opm/autodiff/WellInterface_impl.hpp @@ -979,6 +979,48 @@ namespace Opm + + template + bool + WellInterface::isVFPActive() const + { + // since the well_controls only handles the VFP number when THP constraint/target is there. + // we need to get the table number through the parser, in case THP constraint/target is not there. + // When THP control/limit is not active, if available VFP table is provided, we will still need to + // update THP value. However, it will only used for output purpose. + + if (well_type_ == PRODUCER) { // producer + const int table_id = well_ecl_->getProductionProperties(current_step_).VFPTableNumber; + if (table_id <= 0) { + return false; + } else { + if (vfp_properties_->getProd()->hasTable(table_id)) { + return true; + } else { + OPM_THROW(std::runtime_error, "VFPPROD table " << std::to_string(table_id) << " is specfied," + << " for well " << name() << ", while we could not access it during simulation"); + } + } + + } else { // injector + const int table_id = well_ecl_->getInjectionProperties(current_step_).VFPTableNumber; + if (table_id <= 0) { + return false; + } else { + if (vfp_properties_->getInj()->hasTable(table_id)) { + return true; + } else { + OPM_THROW(std::runtime_error, "VFPINJ table " << std::to_string(table_id) << " is specfied," + << " for well " << name() << ", while we could not access it during simulation"); + } + } + } + } + + + + + template void WellInterface::calculateReservoirRates(WellState& well_state) const