diff --git a/opm/autodiff/BlackoilModelParametersEbos.hpp b/opm/autodiff/BlackoilModelParametersEbos.hpp index f9a86bd30..b42a2aa1e 100644 --- a/opm/autodiff/BlackoilModelParametersEbos.hpp +++ b/opm/autodiff/BlackoilModelParametersEbos.hpp @@ -72,7 +72,7 @@ SET_BOOL_PROP(FlowModelParameters, UpdateEquationsScaling, false); SET_BOOL_PROP(FlowModelParameters, UseUpdateStabilization, true); SET_BOOL_PROP(FlowModelParameters, MatrixAddWellContributions, false); SET_SCALAR_PROP(FlowModelParameters, TolerancePressureMsWells, 0.01 *1e5); -SET_SCALAR_PROP(FlowModelParameters, MaxPressureChangeMsWells, 2.0 *1e5); +SET_SCALAR_PROP(FlowModelParameters, MaxPressureChangeMsWells, 1e6); SET_BOOL_PROP(FlowModelParameters, UseInnerIterationsMsWells, true); SET_INT_PROP(FlowModelParameters, MaxInnerIterMsWells, 100); diff --git a/opm/simulators/wells/BlackoilWellModel.hpp b/opm/simulators/wells/BlackoilWellModel.hpp index 977a61354..3e74cfe49 100644 --- a/opm/simulators/wells/BlackoilWellModel.hpp +++ b/opm/simulators/wells/BlackoilWellModel.hpp @@ -281,7 +281,7 @@ namespace Opm { std::vector is_cell_perforated_; // create the well container - std::vector createWellContainer(const int time_step, Opm::DeferredLogger& deferred_logger); + std::vector createWellContainer(const int time_step, const Wells* wells, const bool allow_closing_opening_wells, Opm::DeferredLogger& deferred_logger); WellInterfacePtr createWellForWellTest(const std::string& well_name, const int report_step, Opm::DeferredLogger& deferred_logger) const; @@ -337,9 +337,11 @@ namespace Opm { // xw to update Well State void recoverWellSolutionAndUpdateWellState(const BVector& x); - void updateWellControls(Opm::DeferredLogger& deferred_logger); + void updateWellControls(const std::vector& B_avg, + Opm::DeferredLogger& deferred_logger); - void updateGroupControls(Opm::DeferredLogger& deferred_logger); + void updateGroupControls(const std::vector& B_avg, + Opm::DeferredLogger& deferred_logger); // setting the well_solutions_ based on well_state. void updatePrimaryVariables(Opm::DeferredLogger& deferred_logger); @@ -387,7 +389,8 @@ namespace Opm { // some preparation work, mostly related to group control and RESV, // at the beginning of each time step (Not report step) - void prepareTimeStep(Opm::DeferredLogger& deferred_logger); + void prepareTimeStep(const std::vector& B_avg, + Opm::DeferredLogger& deferred_logger); void prepareGroupControl(Opm::DeferredLogger& deferred_logger); diff --git a/opm/simulators/wells/BlackoilWellModel_impl.hpp b/opm/simulators/wells/BlackoilWellModel_impl.hpp index 64ed7b096..b9c076d06 100644 --- a/opm/simulators/wells/BlackoilWellModel_impl.hpp +++ b/opm/simulators/wells/BlackoilWellModel_impl.hpp @@ -306,7 +306,7 @@ namespace Opm { wellTesting(reportStepIdx, simulationTime, local_deferredLogger); // create the well container - well_container_ = createWellContainer(reportStepIdx, local_deferredLogger); + well_container_ = createWellContainer(reportStepIdx, wells(), /*allow_closing_opening_wells=*/true, local_deferredLogger); // do the initialization for all the wells // TODO: to see whether we can postpone of the intialization of the well containers to @@ -414,7 +414,6 @@ namespace Opm { updateWellTestState(simulationTime, wellTestState_); // calculate the well potentials for output - // TODO: when necessary try { std::vector well_potentials; computeWellPotentials(well_potentials, local_deferredLogger); @@ -513,7 +512,7 @@ namespace Opm { template std::vector::WellInterfacePtr > BlackoilWellModel:: - createWellContainer(const int time_step, Opm::DeferredLogger& deferred_logger) + createWellContainer(const int time_step, const Wells* wells, const bool allow_closing_opening_wells, Opm::DeferredLogger& deferred_logger) { std::vector well_container; @@ -525,7 +524,7 @@ namespace Opm { // With the following way, it will have the same order with wells struct // Hopefully, it can generate the same residual history with master branch for (int w = 0; w < nw; ++w) { - const std::string well_name = std::string(wells()->name[w]); + const std::string well_name = std::string(wells->name[w]); // finding the location of the well in wells_ecl const int nw_wells_ecl = wells_ecl_.size(); @@ -543,58 +542,61 @@ namespace Opm { const Well2& well_ecl = wells_ecl_[index_well]; - // A new WCON keywords can re-open a well that was closed/shut due to Physical limit - if ( wellTestState_.hasWell(well_name, WellTestConfig::Reason::PHYSICAL ) ) { - // TODO: more checking here, to make sure this standard more specific and complete - // maybe there is some WCON keywords will not open the well - if (well_state_.effectiveEventsOccurred(w)) { - if (wellTestState_.lastTestTime(well_name) == ebosSimulator_.time()) { - // The well was shut this timestep, we are most likely retrying - // a timestep without the well in question, after it caused - // repeated timestep cuts. It should therefore not be opened, - // even if it was new or received new targets this report step. - well_state_.setEffectiveEventsOccurred(w, false); - } else { - wellTestState_.openWell(well_name); + if (allow_closing_opening_wells) { + // A new WCON keywords can re-open a well that was closed/shut due to Physical limit + if ( wellTestState_.hasWell(well_name, WellTestConfig::Reason::PHYSICAL ) ) { + // TODO: more checking here, to make sure this standard more specific and complete + // maybe there is some WCON keywords will not open the well + if (well_state_.effectiveEventsOccurred(w)) { + if (wellTestState_.lastTestTime(well_name) == ebosSimulator_.time()) { + // The well was shut this timestep, we are most likely retrying + // a timestep without the well in question, after it caused + // repeated timestep cuts. It should therefore not be opened, + // even if it was new or received new targets this report step. + well_state_.setEffectiveEventsOccurred(w, false); + } else { + wellTestState_.openWell(well_name); + } } } - } - // TODO: should we do this for all kinds of closing reasons? - // something like wellTestState_.hasWell(well_name)? - if ( wellTestState_.hasWell(well_name, WellTestConfig::Reason::ECONOMIC) || - wellTestState_.hasWell(well_name, WellTestConfig::Reason::PHYSICAL) ) { - if( well_ecl.getAutomaticShutIn() ) { - // shut wells are not added to the well container - // TODO: make a function from well_state side to handle the following - well_state_.thp()[w] = 0.; - well_state_.bhp()[w] = 0.; - const int np = numPhases(); - for (int p = 0; p < np; ++p) { - well_state_.wellRates()[np * w + p] = 0.; + + // TODO: should we do this for all kinds of closing reasons? + // something like wellTestState_.hasWell(well_name)? + if ( wellTestState_.hasWell(well_name, WellTestConfig::Reason::ECONOMIC) || + wellTestState_.hasWell(well_name, WellTestConfig::Reason::PHYSICAL) ) { + if( well_ecl.getAutomaticShutIn() ) { + // shut wells are not added to the well container + // TODO: make a function from well_state side to handle the following + well_state_.thp()[w] = 0.; + well_state_.bhp()[w] = 0.; + const int np = numPhases(); + for (int p = 0; p < np; ++p) { + well_state_.wellRates()[np * w + p] = 0.; + } + continue; + } else { + // close wells are added to the container but marked as closed + struct WellControls* well_controls = wells->ctrls[w]; + well_controls_stop_well(well_controls); } - continue; - } else { - // close wells are added to the container but marked as closed - struct WellControls* well_controls = wells()->ctrls[w]; - well_controls_stop_well(well_controls); } } // Use the pvtRegionIdx from the top cell - const int well_cell_top = wells()->well_cells[wells()->well_connpos[w]]; + const int well_cell_top = wells->well_cells[wells->well_connpos[w]]; const int pvtreg = pvt_region_idx_[well_cell_top]; if ( !well_ecl.isMultiSegment() || !param_.use_multisegment_well_) { if ( GET_PROP_VALUE(TypeTag, EnablePolymerMW) && well_ecl.isInjector() ) { - well_container.emplace_back(new StandardWellV(well_ecl, time_step, wells(), + well_container.emplace_back(new StandardWellV(well_ecl, time_step, wells, param_, *rateConverter_, pvtreg, numComponents() ) ); } else { - well_container.emplace_back(new StandardWell(well_ecl, time_step, wells(), + well_container.emplace_back(new StandardWell(well_ecl, time_step, wells, param_, *rateConverter_, pvtreg, numComponents() ) ); } } else { - well_container.emplace_back(new MultisegmentWell(well_ecl, time_step, wells(), + well_container.emplace_back(new MultisegmentWell(well_ecl, time_step, wells, param_, *rateConverter_, pvtreg, numComponents() ) ); } } @@ -679,18 +681,18 @@ namespace Opm { int exception_thrown = 0; try { - if (iterationIdx == 0) { - calculateExplicitQuantities(local_deferredLogger); - prepareTimeStep(local_deferredLogger); - } - - updateWellControls(local_deferredLogger); - // Set the well primary variables based on the value of well solutions - initPrimaryVariablesEvaluation(); - std::vector< Scalar > B_avg(numComponents(), Scalar() ); computeAverageFormationFactor(B_avg); + if (iterationIdx == 0) { + calculateExplicitQuantities(local_deferredLogger); + prepareTimeStep(B_avg, local_deferredLogger); + } + + updateWellControls(B_avg, local_deferredLogger); + // Set the well primary variables based on the value of well solutions + initPrimaryVariablesEvaluation(); + if (param_.solve_welleq_initially_ && iterationIdx == 0) { // solve the well equations as a pre-processing step last_report_ = solveWellEq(B_avg, dt, local_deferredLogger); @@ -699,7 +701,7 @@ namespace Opm { if (initial_step_) { // update the explicit quantities to get the initial fluid distribution in the well correct. calculateExplicitQuantities(local_deferredLogger); - prepareTimeStep(local_deferredLogger); + prepareTimeStep(B_avg, local_deferredLogger); last_report_ = solveWellEq(B_avg, dt, local_deferredLogger); initial_step_ = false; } @@ -708,7 +710,7 @@ namespace Opm { // reservoir state, will tihs be a better place to inialize the explict information? } - assembleWellEq(B_avg, dt, local_deferredLogger); + assembleWellEq(B_avg, dt, local_deferredLogger); } catch (std::exception& e) { exception_thrown = 1; @@ -921,7 +923,7 @@ namespace Opm { // are active wells anywhere in the global domain. if( wellsActive() ) { - updateWellControls(deferred_logger); + updateWellControls(B_avg, deferred_logger); initPrimaryVariablesEvaluation(); } } catch (std::exception& e) { @@ -1025,7 +1027,8 @@ namespace Opm { template void BlackoilWellModel:: - updateWellControls(Opm::DeferredLogger& deferred_logger) + updateWellControls(const std::vector& B_avg, + Opm::DeferredLogger& deferred_logger) { // Even if there are no wells active locally, we cannot // return as the DeferredLogger uses global communication. @@ -1033,10 +1036,10 @@ namespace Opm { if( !wellsActive() ) return ; for (const auto& well : well_container_) { - well->updateWellControl(ebosSimulator_, well_state_, deferred_logger); + well->updateWellControl(ebosSimulator_, B_avg, well_state_, deferred_logger); } - updateGroupControls(deferred_logger); + updateGroupControls(B_avg, deferred_logger); } @@ -1066,17 +1069,74 @@ namespace Opm { BlackoilWellModel:: computeWellPotentials(std::vector& well_potentials, Opm::DeferredLogger& deferred_logger) { - Opm::DeferredLogger local_deferredLogger; // number of wells and phases const int nw = numWells(); const int np = numPhases(); well_potentials.resize(nw * np, 0.0); + + const int reportStepIdx = ebosSimulator_.episodeIndex(); + const double invalid_alq = -1e100; + const double invalid_vfp = -2147483647; + auto well_state_copy = well_state_; + const Wells* local_wells = clone_wells(wells()); + std::vector well_container_copy = createWellContainer(reportStepIdx, local_wells, /*allow_closing_opening_wells=*/ false, deferred_logger); + + // average B factors are required for the convergence checking of well equations + // Note: this must be done on all processes, even those with + // no wells needing testing, otherwise we will have locking. + std::vector< Scalar > B_avg(numComponents(), Scalar() ); + computeAverageFormationFactor(B_avg); + const Opm::SummaryConfig& summaryConfig = ebosSimulator_.vanguard().summaryConfig(); int exception_thrown = 0; try { - for (const auto& well : well_container_) { + for (const auto& well : well_container_copy) { // Only compute the well potential when asked for + well->init(&phase_usage_, depth_, gravity_, number_of_cells_); + + WellControls* wc = well->wellControls(); + well_controls_clear(wc); + well_controls_assert_number_of_phases( wc , np); + if (well->wellType() == INJECTOR) { + const auto& injectionProperties = well->wellEcl()->getInjectionProperties(); + + if (injectionProperties.hasInjectionControl(WellInjector::THP)) { + const double thp_limit = injectionProperties.THPLimit; + const int vfp_number = injectionProperties.VFPTableNumber; + well_controls_add_new(THP, thp_limit, invalid_alq, vfp_number, NULL, wc); + } + + // we always have a bhp limit + const double bhp_limit = injectionProperties.BHPLimit; + well_controls_add_new(BHP, bhp_limit, invalid_alq, invalid_vfp, NULL, wc); + } else { + const auto& productionProperties = well->wellEcl()->getProductionProperties(); + if (productionProperties.hasProductionControl(WellProducer::THP)) { + const double thp_limit = productionProperties.THPLimit; + const double alq_value = productionProperties.ALQValue; + const int vfp_number = productionProperties.VFPTableNumber; + well_controls_add_new(THP, thp_limit, alq_value, vfp_number, NULL, wc); + } + + // we always have a bhp limit + const double bhp_limit = productionProperties.BHPLimit; + well_controls_add_new(BHP, bhp_limit, invalid_alq, invalid_vfp, NULL, wc); + + well->setVFPProperties(vfp_properties_.get()); + + } + + if (has_polymer_) + { + const Grid& grid = ebosSimulator_.vanguard().grid(); + if (PolymerModule::hasPlyshlog() || GET_PROP_VALUE(TypeTag, EnablePolymerMW) ) { + well->computeRepRadiusPerfLength(grid, cartesian_to_compressed_, deferred_logger); + } + } + + + const bool needed_for_output = ((summaryConfig.hasSummaryKey( "WWPI:" + well->name()) || summaryConfig.hasSummaryKey( "WOPI:" + well->name()) || summaryConfig.hasSummaryKey( "WGPI:" + well->name())) && well->wellType() == INJECTOR) || @@ -1086,7 +1146,7 @@ namespace Opm { if (needed_for_output || wellCollection().requireWellPotentials()) { std::vector potentials; - well->computeWellPotentials(ebosSimulator_, well_state_, potentials, deferred_logger); + well->computeWellPotentials(ebosSimulator_, B_avg, well_state_copy, potentials, deferred_logger); // putting the sucessfully calculated potentials to the well_potentials for (int p = 0; p < np; ++p) { well_potentials[well->indexOfWell() * np + p] = std::abs(potentials[p]); @@ -1111,7 +1171,8 @@ namespace Opm { template void BlackoilWellModel:: - prepareTimeStep(Opm::DeferredLogger& deferred_logger) + prepareTimeStep(const std::vector& B_avg, + Opm::DeferredLogger& deferred_logger) { if ( wellCollection().havingVREPGroups() ) { @@ -1132,7 +1193,7 @@ namespace Opm { int exception_thrown = 0; try { for (const auto& well : well_container_) { - well->checkWellOperability(ebosSimulator_, well_state_, deferred_logger); + well->checkWellOperability(ebosSimulator_, B_avg, well_state_, deferred_logger); } // since the controls are all updated, we should update well_state accordingly for (const auto& well : well_container_) { @@ -1144,7 +1205,7 @@ namespace Opm { if (!well->isOperable() ) continue; if (well_state_.effectiveEventsOccurred(w) ) { - well->updateWellStateWithTarget(ebosSimulator_, well_state_, deferred_logger); + well->updateWellStateWithTarget(ebosSimulator_, B_avg, well_state_, deferred_logger); } // there is no new well control change input within a report step, @@ -1372,7 +1433,8 @@ namespace Opm { template void BlackoilWellModel:: - updateGroupControls(Opm::DeferredLogger& deferred_logger) + updateGroupControls(const std::vector& B_avg, + Opm::DeferredLogger& deferred_logger) { if (wellCollection().groupControlActive()) { @@ -1399,7 +1461,7 @@ namespace Opm { // TODO: we should only do the well is involved in the update group targets for (auto& well : well_container_) { - well->updateWellStateWithTarget(ebosSimulator_, well_state_, deferred_logger); + well->updateWellStateWithTarget(ebosSimulator_, B_avg, well_state_, deferred_logger); well->updatePrimaryVariables(well_state_, deferred_logger); } } diff --git a/opm/simulators/wells/MultisegmentWell.hpp b/opm/simulators/wells/MultisegmentWell.hpp index b7eaae8e1..aeeee9050 100644 --- a/opm/simulators/wells/MultisegmentWell.hpp +++ b/opm/simulators/wells/MultisegmentWell.hpp @@ -120,8 +120,9 @@ namespace Opm /// updating the well state based the current control mode virtual void updateWellStateWithTarget(const Simulator& ebos_simulator, + const std::vector& B_avg, WellState& well_state, - Opm::DeferredLogger& deferred_logger) const override; + Opm::DeferredLogger& deferred_logger) override; /// check whether the well equations get converged for this well virtual ConvergenceReport getWellConvergence(const std::vector& B_avg, Opm::DeferredLogger& deferred_logger) const override; @@ -139,6 +140,7 @@ namespace Opm /// computing the well potentials for group control virtual void computeWellPotentials(const Simulator& ebosSimulator, + const std::vector& B_avg, const WellState& well_state, std::vector& well_potentials, Opm::DeferredLogger& deferred_logger) override; @@ -305,6 +307,8 @@ namespace Opm const bool& allow_cf, std::vector& cq_s, EvalWell& perf_press, + double& perf_dis_gas_rate, + double& perf_vap_oil_rate, Opm::DeferredLogger& deferred_logger) const; // convert a Eval from reservoir to contain the derivative related to wells @@ -327,6 +331,13 @@ namespace Opm const int perf, std::vector& mob) const; + virtual void computeWellRatesWithBhp(const Simulator& ebosSimulator, + const std::vector& B_avg, + const double& bhp, + const bool iterate, + std::vector& well_flux, + Opm::DeferredLogger& deferred_logger) override; + void assembleControlEq(Opm::DeferredLogger& deferred_logger) const; void assemblePressureEq(const int seg) const; @@ -345,6 +356,7 @@ namespace Opm // checking the operability of the well based on current reservoir condition // it is not implemented for multisegment well yet virtual void checkWellOperability(const Simulator& ebos_simulator, + const std::vector& B_avg, const WellState& well_state, Opm::DeferredLogger& deferred_logger) override; @@ -366,7 +378,7 @@ namespace Opm WellState& well_state, Opm::DeferredLogger& deferred_logger); - virtual void wellTestingPhysical(Simulator& simulator, const std::vector& B_avg, + virtual void wellTestingPhysical(const Simulator& simulator, const std::vector& B_avg, const double simulation_time, const int report_step, WellState& well_state, WellTestState& welltest_state, Opm::DeferredLogger& deferred_logger) override; diff --git a/opm/simulators/wells/MultisegmentWell_impl.hpp b/opm/simulators/wells/MultisegmentWell_impl.hpp index bd6328a0f..580c50483 100644 --- a/opm/simulators/wells/MultisegmentWell_impl.hpp +++ b/opm/simulators/wells/MultisegmentWell_impl.hpp @@ -260,8 +260,9 @@ namespace Opm void MultisegmentWell:: updateWellStateWithTarget(const Simulator& /* ebos_simulator */, + const std::vector& B_avg, WellState& well_state, - Opm::DeferredLogger& /* deferred_logger */) const + Opm::DeferredLogger& /* deferred_logger */) { // Updating well state bas on well control // Target values are used as initial conditions for BHP, THP, and SURFACE_RATE @@ -541,24 +542,73 @@ namespace Opm template void MultisegmentWell:: - computeWellPotentials(const Simulator& /* ebosSimulator */, - const WellState& /* well_state */, + computeWellPotentials(const Simulator& ebosSimulator, + const std::vector& B_avg, + const WellState& well_state, std::vector& well_potentials, Opm::DeferredLogger& deferred_logger) { - const std::string msg = std::string("Well potential calculation is not supported for multisegment wells \n") - + "A well potential of zero is returned for output purposes. \n" - + "If you need well potential to set the guide rate for group controled wells \n" - + "you will have to change the " + name() + " well to a standard well \n"; - - deferred_logger.warning("WELL_POTENTIAL_NOT_IMPLEMENTED_FOR_MULTISEG_WELLS", msg); - const int np = number_of_phases_; well_potentials.resize(np, 0.0); + updatePrimaryVariables(well_state, deferred_logger); + + // initialize the primary variables in Evaluation, which is used in computePerfRate for computeWellPotentials + // TODO: for computeWellPotentials, no derivative is required actually + initPrimaryVariablesEvaluation(); + + // get the bhp value based on the bhp constraints + const double bhp = Base::mostStrictBhpFromBhpLimits(deferred_logger); + + // does the well have a THP related constraint? + if ( !Base::wellHasTHPConstraints() ) { + assert(std::abs(bhp) != std::numeric_limits::max()); + + computeWellRatesWithBhp(ebosSimulator, B_avg, bhp, /*iterate=*/ true, well_potentials, deferred_logger); + } else { + + const std::string msg = std::string("Well potential calculation is not supported for thp controlled multisegment wells \n") + + "A well potential of zero is returned for output purposes. \n" + + "If you need well potential computed from thp to set the guide rate for group controled wells \n" + + "you will have to change the " + name() + " well to a standard well \n"; + + deferred_logger.warning("WELL_POTENTIAL_FOR_THP_NOT_IMPLEMENTED_FOR_MULTISEG_WELLS", msg); + return; + } + } + template + void + MultisegmentWell:: + computeWellRatesWithBhp(const Simulator& ebosSimulator, + const std::vector& B_avg, + const double& bhp, + const bool iterate, + std::vector& well_flux, + Opm::DeferredLogger& deferred_logger) + { + // store a copy of the well state, we don't want to update the real well state + WellState copy = ebosSimulator.problem().wellModel().wellState(); + + initPrimaryVariablesEvaluation(); + + if (iterate) { + const double dt = ebosSimulator.timeStepSize(); + // iterate to get a solution that satisfies the bhp potential. + iterateWellEquations(ebosSimulator, B_avg, dt, copy, deferred_logger); + } + + // compute the potential and store in the flux vector. + const int np = number_of_phases_; + well_flux.resize(np, 0.0); + for(int compIdx = 0; compIdx < num_components_; ++compIdx) { + const EvalWell rate = getSegmentRate(0, compIdx); + well_flux[ebosCompIdxToFlowCompIdx(compIdx)] += rate.value(); + } + + } @@ -743,7 +793,7 @@ namespace Opm MultisegmentWell:: updateWellState(const BVectorWell& dwells, WellState& well_state, - Opm::DeferredLogger& /* deferred_logger */, + Opm::DeferredLogger& deferred_logger, const double relaxation_factor) const { const double dFLimit = param_.dwell_fraction_max_; @@ -773,7 +823,7 @@ namespace Opm const int sign = dwells[seg][SPres] > 0.? 1 : -1; const double dx_limited = sign * std::min(std::abs(dwells[seg][SPres]), relaxation_factor * max_pressure_change); // const double dx_limited = sign * std::min(std::abs(dwells[seg][SPres]) * relaxation_factor, max_pressure_change); - primary_variables_[seg][SPres] = old_primary_variables[seg][SPres] - dx_limited; + primary_variables_[seg][SPres] = std::max( old_primary_variables[seg][SPres] - dx_limited, 1e5); } // update the total rate // TODO: should we have a limitation of the total rate change? @@ -961,6 +1011,8 @@ namespace Opm const bool& allow_cf, std::vector& cq_s, EvalWell& perf_press, + double& perf_dis_gas_rate, + double& perf_vap_oil_rate, Opm::DeferredLogger& deferred_logger) const { @@ -1077,6 +1129,29 @@ namespace Opm cq_s[comp_idx] = cmix_s[comp_idx] * cqt_is; } } // end for injection perforations + + // calculating the perforation solution gas rate and solution oil rates + if (well_type_ == PRODUCER) { + if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) { + const unsigned oilCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx); + const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx); + // TODO: the formulations here remain to be tested with cases with strong crossflow through production wells + // s means standard condition, r means reservoir condition + // q_os = q_or * b_o + rv * q_gr * b_g + // q_gs = q_gr * g_g + rs * q_or * b_o + // d = 1.0 - rs * rv + // q_or = 1 / (b_o * d) * (q_os - rv * q_gs) + // q_gr = 1 / (b_g * d) * (q_gs - rs * q_os) + + const double d = 1.0 - rv.value() * rs.value(); + // vaporized oil into gas + // rv * q_gr * b_g = rv * (q_gs - rs * q_os) / d + perf_vap_oil_rate = rv.value() * (cq_s[gasCompIdx].value() - rs.value() * cq_s[oilCompIdx].value()) / d; + // dissolved of gas in oil + // rs * q_or * b_o = rs * (q_os - rv * q_gs) / d + perf_dis_gas_rate = rs.value() * (cq_s[oilCompIdx].value() - rv.value() * cq_s[gasCompIdx].value()) / d; + } + } } @@ -1673,6 +1748,7 @@ namespace Opm void MultisegmentWell:: checkWellOperability(const Simulator& /* ebos_simulator */, + const std::vector& /*B_avg */, const WellState& /* well_state */, Opm::DeferredLogger& deferred_logger) { @@ -1829,7 +1905,7 @@ namespace Opm updateWellState(dx_well, well_state, deferred_logger, relaxation_factor); // TODO: should we do something more if a switching of control happens - this->updateWellControl(ebosSimulator, well_state, deferred_logger); + this->updateWellControl(ebosSimulator, B_avg, well_state, deferred_logger); initPrimaryVariablesEvaluation(); } @@ -1877,6 +1953,9 @@ namespace Opm duneD_ = 0.0; resWell_ = 0.0; + well_state.wellVaporizedOilRates()[index_of_well_] = 0.; + well_state.wellDissolvedGasRates()[index_of_well_] = 0.; + // for the black oil cases, there will be four equations, // the first three of them are the mass balance equations, the last one is the pressure equations. // @@ -1944,7 +2023,16 @@ namespace Opm getMobility(ebosSimulator, perf, mob); std::vector cq_s(num_components_, 0.0); EvalWell perf_press; - computePerfRatePressure(int_quants, mob, seg, perf, seg_pressure, allow_cf, cq_s, perf_press, deferred_logger); + double perf_dis_gas_rate = 0.; + double perf_vap_oil_rate = 0.; + + computePerfRatePressure(int_quants, mob, seg, perf, seg_pressure, allow_cf, cq_s, perf_press, perf_dis_gas_rate, perf_vap_oil_rate, deferred_logger); + + // updating the solution gas rate and solution oil rate + if (well_type_ == PRODUCER) { + well_state.wellDissolvedGasRates()[index_of_well_] += perf_dis_gas_rate; + well_state.wellVaporizedOilRates()[index_of_well_] += perf_vap_oil_rate; + } // store the perf pressure and rates const int rate_start_offset = (first_perf_ + perf) * number_of_phases_; @@ -1995,7 +2083,7 @@ namespace Opm template void MultisegmentWell:: - wellTestingPhysical(Simulator& /* simulator */, const std::vector& /* B_avg */, + wellTestingPhysical(const Simulator& /* simulator */, const std::vector& /* B_avg */, const double /* simulation_time */, const int /* report_step */, WellState& /* well_state */, WellTestState& /* welltest_state */, Opm::DeferredLogger& deferred_logger) { diff --git a/opm/simulators/wells/StandardWell.hpp b/opm/simulators/wells/StandardWell.hpp index e52572f37..bc68d054e 100644 --- a/opm/simulators/wells/StandardWell.hpp +++ b/opm/simulators/wells/StandardWell.hpp @@ -148,8 +148,9 @@ namespace Opm Opm::DeferredLogger& deferred_logger) override; virtual void updateWellStateWithTarget(const Simulator& ebos_simulator, + const std::vector& B_avg, WellState& well_state, - Opm::DeferredLogger& deferred_logger) const override; + Opm::DeferredLogger& deferred_logger) override; /// check whether the well equations get converged for this well virtual ConvergenceReport getWellConvergence(const std::vector& B_avg, @@ -168,6 +169,7 @@ namespace Opm /// computing the well potentials for group control virtual void computeWellPotentials(const Simulator& ebosSimulator, + const std::vector& B_avg, const WellState& well_state, std::vector& well_potentials, Opm::DeferredLogger& deferred_logger) /* const */ override; @@ -327,17 +329,18 @@ namespace Opm double& perf_vap_oil_rate, Opm::DeferredLogger& deferred_logger) const; - // TODO: maybe we should provide a light version of computePerfRate, which does not include the - // calculation of the derivatives - void computeWellRatesWithBhp(const Simulator& ebosSimulator, - const EvalWell& bhp, - std::vector& well_flux, - Opm::DeferredLogger& deferred_logger) const; + virtual void computeWellRatesWithBhp(const Simulator& ebosSimulator, + const std::vector& B_avg, + const double& bhp, + const bool iterate, + std::vector& well_flux, + Opm::DeferredLogger& deferred_logger) override; std::vector computeWellPotentialWithTHP(const Simulator& ebosSimulator, + const std::vector& B_avg, const double initial_bhp, // bhp from BHP constraints const std::vector& initial_potential, - Opm::DeferredLogger& deferred_logger) const; + Opm::DeferredLogger& deferred_logger); template ValueType calculateBhpFromThp(const std::vector& rates, const int control_index, Opm::DeferredLogger& deferred_logger) const; @@ -373,6 +376,7 @@ namespace Opm // update the operability status of the well is operable under the current reservoir condition // mostly related to BHP limit and THP limit virtual void checkWellOperability(const Simulator& ebos_simulator, + const std::vector& B_avg, const WellState& well_state, Opm::DeferredLogger& deferred_logger ) override; @@ -380,24 +384,29 @@ namespace Opm // check whether the well is operable under the current reservoir condition // mostly related to BHP limit and THP limit void updateWellOperability(const Simulator& ebos_simulator, + const std::vector& B_avg, const WellState& well_state, Opm::DeferredLogger& deferred_logger ); // check whether the well is operable under BHP limit with current reservoir condition - void checkOperabilityUnderBHPLimitProducer(const Simulator& ebos_simulator, Opm::DeferredLogger& deferred_logger); + void checkOperabilityUnderBHPLimitProducer(const Simulator& ebos_simulator, + const std::vector& B_avg, + Opm::DeferredLogger& deferred_logger); // check whether the well is operable under THP limit with current reservoir condition void checkOperabilityUnderTHPLimitProducer(const Simulator& ebos_simulator, Opm::DeferredLogger& deferred_logger); // update WellState based on IPR and associated VFP table void updateWellStateWithTHPTargetIPR(const Simulator& ebos_simulator, + const std::vector& B_avg, WellState& well_state, - Opm::DeferredLogger& deferred_logger) const; + Opm::DeferredLogger& deferred_logger); void updateWellStateWithTHPTargetIPRProducer(const Simulator& ebos_simulator, + const std::vector& B_avg, WellState& well_state, - Opm::DeferredLogger& deferred_logger) const; + Opm::DeferredLogger& deferred_logger); // for a well, when all drawdown are in the wrong direction, then this well will not // be able to produce/inject . @@ -433,7 +442,7 @@ namespace Opm static double relaxationFactorRate(const std::vector& primary_variables, const BVectorWell& dwells); - virtual void wellTestingPhysical(Simulator& simulator, const std::vector& B_avg, + virtual void wellTestingPhysical(const Simulator& simulator, const std::vector& B_avg, const double simulation_time, const int report_step, WellState& well_state, WellTestState& welltest_state, Opm::DeferredLogger& deferred_logger) override; diff --git a/opm/simulators/wells/StandardWellV.hpp b/opm/simulators/wells/StandardWellV.hpp index 1cbdfd129..246244702 100644 --- a/opm/simulators/wells/StandardWellV.hpp +++ b/opm/simulators/wells/StandardWellV.hpp @@ -149,8 +149,9 @@ namespace Opm Opm::DeferredLogger& deferred_logger) override; virtual void updateWellStateWithTarget(const Simulator& ebos_simulator, + const std::vector& B_avg, WellState& well_state, - Opm::DeferredLogger& deferred_logger) const override; + Opm::DeferredLogger& deferred_logger) override; /// check whether the well equations get converged for this well virtual ConvergenceReport getWellConvergence(const std::vector& B_avg, Opm::DeferredLogger& deferred_logger) const override; @@ -168,6 +169,7 @@ namespace Opm /// computing the well potentials for group control virtual void computeWellPotentials(const Simulator& ebosSimulator, + const std::vector& B_avg, const WellState& well_state, std::vector& well_potentials, Opm::DeferredLogger& deferred_logger) /* const */ override; @@ -330,17 +332,18 @@ namespace Opm double& perf_vap_oil_rate, Opm::DeferredLogger& deferred_logger) const; - // TODO: maybe we should provide a light version of computePerfRate, which does not include the - // calculation of the derivatives - void computeWellRatesWithBhp(const Simulator& ebosSimulator, - const EvalWell& bhp, - std::vector& well_flux, - Opm::DeferredLogger& deferred_logger) const; + virtual void computeWellRatesWithBhp(const Simulator& ebosSimulator, + const std::vector& B_avg, + const double& bhp, + const bool iterate, + std::vector& well_flux, + Opm::DeferredLogger& deferred_logger) override; std::vector computeWellPotentialWithTHP(const Simulator& ebosSimulator, + const std::vector& B_avg, const double initial_bhp, // bhp from BHP constraints const std::vector& initial_potential, - Opm::DeferredLogger& deferred_logger) const; + Opm::DeferredLogger& deferred_logger); template ValueType calculateBhpFromThp(const std::vector& rates, const int control_index, Opm::DeferredLogger& deferred_logger) const; @@ -376,29 +379,35 @@ namespace Opm // update the operability status of the well is operable under the current reservoir condition // mostly related to BHP limit and THP limit virtual void checkWellOperability(const Simulator& ebos_simulator, + const std::vector& B_avg, const WellState& well_state, Opm::DeferredLogger& deferred_logger) override; // check whether the well is operable under the current reservoir condition // mostly related to BHP limit and THP limit void updateWellOperability(const Simulator& ebos_simulator, + const std::vector& B_avg, const WellState& well_state, Opm::DeferredLogger& deferred_logger); // check whether the well is operable under BHP limit with current reservoir condition - void checkOperabilityUnderBHPLimitProducer(const Simulator& ebos_simulator, Opm::DeferredLogger& deferred_logger); + void checkOperabilityUnderBHPLimitProducer(const Simulator& ebos_simulator, + const std::vector& B_avg, + Opm::DeferredLogger& deferred_logger); // check whether the well is operable under THP limit with current reservoir condition void checkOperabilityUnderTHPLimitProducer(const Simulator& ebos_simulator, Opm::DeferredLogger& deferred_logger); // update WellState based on IPR and associated VFP table void updateWellStateWithTHPTargetIPR(const Simulator& ebos_simulator, + const std::vector& B_avg, WellState& well_state, - Opm::DeferredLogger& deferred_logger) const; + Opm::DeferredLogger& deferred_logger); void updateWellStateWithTHPTargetIPRProducer(const Simulator& ebos_simulator, + const std::vector& B_avg, WellState& well_state, - Opm::DeferredLogger& deferred_logger) const; + Opm::DeferredLogger& deferred_logger); // for a well, when all drawdown are in the wrong direction, then this well will not // be able to produce/inject . @@ -434,7 +443,7 @@ namespace Opm static double relaxationFactorRate(const std::vector& primary_variables, const BVectorWell& dwells); - virtual void wellTestingPhysical(Simulator& simulator, const std::vector& B_avg, + virtual void wellTestingPhysical(const Simulator& simulator, const std::vector& B_avg, const double simulation_time, const int report_step, WellState& well_state, WellTestState& welltest_state, Opm::DeferredLogger& deferred_logger) override; diff --git a/opm/simulators/wells/StandardWellV_impl.hpp b/opm/simulators/wells/StandardWellV_impl.hpp index 3752e2177..b61fe998d 100644 --- a/opm/simulators/wells/StandardWellV_impl.hpp +++ b/opm/simulators/wells/StandardWellV_impl.hpp @@ -480,7 +480,7 @@ namespace Opm void StandardWellV:: assembleWellEq(const Simulator& ebosSimulator, - const std::vector& /* B_avg */, + const std::vector& B_avg, const double dt, WellState& well_state, Opm::DeferredLogger& deferred_logger @@ -488,7 +488,7 @@ namespace Opm { // TODO: only_wells should be put back to save some computation - checkWellOperability(ebosSimulator, well_state, deferred_logger); + checkWellOperability(ebosSimulator, B_avg, well_state, deferred_logger); if (!this->isOperable()) return; @@ -1194,8 +1194,9 @@ namespace Opm void StandardWellV:: updateWellStateWithTarget(const Simulator& ebos_simulator, + const std::vector& B_avg, WellState& well_state, - Opm::DeferredLogger& deferred_logger) const + Opm::DeferredLogger& deferred_logger) { // number of phases const int np = number_of_phases_; @@ -1220,7 +1221,7 @@ namespace Opm case THP: { // when a well can not work under THP target, it switches to BHP control if (this->operability_status_.isOperableUnderTHPLimit() ) { - updateWellStateWithTHPTargetIPR(ebos_simulator, well_state, deferred_logger); + updateWellStateWithTHPTargetIPR(ebos_simulator, B_avg, well_state, deferred_logger); } else { // go to BHP limit assert(this->operability_status_.isOperableUnderBHPLimit() ); @@ -1398,6 +1399,7 @@ namespace Opm void StandardWellV:: checkWellOperability(const Simulator& ebos_simulator, + const std::vector& B_avg, const WellState& well_state, Opm::DeferredLogger& deferred_logger) { @@ -1412,7 +1414,7 @@ namespace Opm const bool old_well_operable = this->operability_status_.isOperable(); - updateWellOperability(ebos_simulator, well_state, deferred_logger); + updateWellOperability(ebos_simulator, B_avg, well_state, deferred_logger); const bool well_operable = this->operability_status_.isOperable(); @@ -1431,6 +1433,7 @@ namespace Opm void StandardWellV:: updateWellOperability(const Simulator& ebos_simulator, + const std::vector& B_avg, const WellState& /* well_state */, Opm::DeferredLogger& deferred_logger) { @@ -1439,7 +1442,7 @@ namespace Opm updateIPR(ebos_simulator, deferred_logger); // checking the BHP limit related - checkOperabilityUnderBHPLimitProducer(ebos_simulator, deferred_logger); + checkOperabilityUnderBHPLimitProducer(ebos_simulator, B_avg, deferred_logger); // checking whether the well can operate under the THP constraints. if (this->wellHasTHPConstraints()) { @@ -1454,7 +1457,9 @@ namespace Opm template void StandardWellV:: - checkOperabilityUnderBHPLimitProducer(const Simulator& ebos_simulator, Opm::DeferredLogger& deferred_logger) + checkOperabilityUnderBHPLimitProducer(const Simulator& ebos_simulator, + const std::vector& B_avg, + Opm::DeferredLogger& deferred_logger) { const double bhp_limit = mostStrictBhpFromBhpLimits(deferred_logger); // Crude but works: default is one atmosphere. @@ -1478,7 +1483,7 @@ namespace Opm // option 2: stick with the above IPR curve // we use IPR here std::vector well_rates_bhp_limit; - computeWellRatesWithBhp(ebos_simulator, EvalWell(numWellEq_ + numEq, bhp_limit), well_rates_bhp_limit, deferred_logger); + computeWellRatesWithBhp(ebos_simulator, B_avg, bhp_limit, /*iterate=*/false, well_rates_bhp_limit, deferred_logger); const double thp = calculateThpFromBhp(well_rates_bhp_limit, bhp_limit, deferred_logger); const double thp_limit = this->getTHPConstraint(deferred_logger); @@ -1583,7 +1588,7 @@ namespace Opm { const double bhp = well_state.bhp()[index_of_well_]; std::vector well_rates; - computeWellRatesWithBhp(ebos_simulator, bhp, well_rates, deferred_logger); + computeWellRatesWithBhp(ebos_simulator, bhp, /*iterate=*/ false, well_rates, deferred_logger); const double sign = (well_type_ == PRODUCER) ? -1. : 1.; const double threshold = sign * std::numeric_limits::min(); @@ -1626,11 +1631,13 @@ namespace Opm void StandardWellV:: updateWellStateWithTHPTargetIPR(const Simulator& ebos_simulator, + const std::vector& B_avg, WellState& well_state, - Opm::DeferredLogger& deferred_logger) const + Opm::DeferredLogger& deferred_logger) { if (well_type_ == PRODUCER) { updateWellStateWithTHPTargetIPRProducer(ebos_simulator, + B_avg, well_state, deferred_logger); } @@ -1650,8 +1657,9 @@ namespace Opm void StandardWellV:: updateWellStateWithTHPTargetIPRProducer(const Simulator& ebos_simulator, + const std::vector& B_avg, WellState& well_state, - Opm::DeferredLogger& deferred_logger) const + Opm::DeferredLogger& deferred_logger) { well_state.thp()[index_of_well_] = this->getTHPConstraint(deferred_logger); @@ -1667,7 +1675,7 @@ namespace Opm initPrimaryVariablesEvaluation(); std::vector rates; - computeWellRatesWithBhp(ebos_simulator, EvalWell(numWellEq_ + numEq, bhp), rates, deferred_logger); + computeWellRatesWithBhp(ebos_simulator, B_avg, bhp, /*iterate=*/false, rates, deferred_logger); // TODO: double checke the obtained rates // this is another places we might obtain negative rates @@ -1690,7 +1698,7 @@ namespace Opm calculateBHPWithTHPTargetIPR(Opm::DeferredLogger& deferred_logger) const { const double thp_target = this->getTHPConstraint(deferred_logger); - const double thp_control_index = this->getTHPControlIndex(); + const double thp_control_index = this->getControlIndex(THP); const int thp_table_id = well_controls_iget_vfp(well_controls_, thp_control_index); const double alq = well_controls_iget_alq(well_controls_, thp_control_index); @@ -2263,12 +2271,42 @@ namespace Opm void StandardWellV:: computeWellRatesWithBhp(const Simulator& ebosSimulator, - const EvalWell& bhp, + const std::vector& B_avg, + const double& bhp, + const bool iterate, std::vector& well_flux, - Opm::DeferredLogger& deferred_logger) const + Opm::DeferredLogger& deferred_logger) { + const int np = number_of_phases_; well_flux.resize(np, 0.0); + WellControls* wc = well_controls_; + const int bhp_index = Base::getControlIndex(BHP); + const double orig_bhp = well_controls_iget_target(wc, bhp_index); + const auto orig_current = well_controls_get_current(wc); + + well_controls_iset_target(wc, bhp_index, bhp); + well_controls_set_current(wc, bhp_index); + + // iterate to get a more accurate well density + if (iterate) { + // create a copy of the well_state to use. If the operability checking is sucessful, we use this one + // to replace the original one + WellState well_state_copy = ebosSimulator.problem().wellModel().wellState(); + well_state_copy.currentControls()[index_of_well_] = bhp_index; + + bool converged = this->solveWellEqUntilConverged(ebosSimulator, B_avg, well_state_copy, deferred_logger); + + if (!converged) { + const std::string msg = " well " + name() + " did not get converged during well potential calculations " + "returning zero values for the potential"; + deferred_logger.debug(msg); + return; + } + updatePrimaryVariables(well_state_copy, deferred_logger); + computeWellConnectionPressures(ebosSimulator, well_state_copy); + initPrimaryVariablesEvaluation(); + } const bool allow_cf = getAllowCrossFlow(); @@ -2278,17 +2316,24 @@ namespace Opm // flux for each perforation std::vector mob(num_components_, {numWellEq_ + numEq, 0.}); getMobility(ebosSimulator, perf, mob, deferred_logger); + //double trans_mult = ebosSimulator.problem().template rockCompTransMultiplier(intQuants, cell_idx); + //const double Tw = well_index_[perf] * trans_mult; std::vector cq_s(num_components_, {numWellEq_ + numEq, 0.}); double perf_dis_gas_rate = 0.; double perf_vap_oil_rate = 0.; - computePerfRate(intQuants, mob, bhp, perf, allow_cf, + computePerfRate(intQuants, mob, EvalWell(numWellEq_ + numEq, bhp), perf, allow_cf, cq_s, perf_dis_gas_rate, perf_vap_oil_rate, deferred_logger); for(int p = 0; p < np; ++p) { well_flux[ebosCompIdxToFlowCompIdx(p)] += cq_s[p].value(); } } + + + // reset bhp limit + well_controls_iset_target(wc, bhp_index, orig_bhp); + well_controls_set_current(wc, orig_current); } @@ -2299,9 +2344,10 @@ namespace Opm std::vector StandardWellV:: computeWellPotentialWithTHP(const Simulator& ebosSimulator, + const std::vector& B_avg, const double initial_bhp, // bhp from BHP constraints const std::vector& initial_potential, - Opm::DeferredLogger& deferred_logger) const + Opm::DeferredLogger& deferred_logger) { // TODO: pay attention to the situation that finally the potential is calculated based on the bhp control // TODO: should we consider the bhp constraints during the iterative process? @@ -2365,7 +2411,7 @@ namespace Opm converged = std::abs(old_bhp - bhp) < bhp_tolerance; - computeWellRatesWithBhp(ebosSimulator, EvalWell(numWellEq_ + numEq, bhp), potentials, deferred_logger); + computeWellRatesWithBhp(ebosSimulator, B_avg, bhp, /*iterate=*/ true, potentials, deferred_logger); // checking whether the potentials have valid values for (const double value : potentials) { @@ -2403,6 +2449,7 @@ namespace Opm void StandardWellV:: computeWellPotentials(const Simulator& ebosSimulator, + const std::vector& B_avg, const WellState& well_state, std::vector& well_potentials, Opm::DeferredLogger& deferred_logger) // const @@ -2424,7 +2471,7 @@ namespace Opm if ( !wellHasTHPConstraints() ) { assert(std::abs(bhp) != std::numeric_limits::max()); - computeWellRatesWithBhp(ebosSimulator, EvalWell(numWellEq_ + numEq, bhp), well_potentials, deferred_logger); + computeWellRatesWithBhp(ebosSimulator, B_avg, bhp, /*iterate=*/ true, well_potentials, deferred_logger); } else { // the well has a THP related constraint // checking whether a well is newly added, it only happens at the beginning of the report step @@ -2436,7 +2483,7 @@ namespace Opm } } else { // We need to generate a reasonable rates to start the iteration process - computeWellRatesWithBhp(ebosSimulator, EvalWell(numWellEq_ + numEq, bhp), well_potentials, deferred_logger); + computeWellRatesWithBhp(ebosSimulator, B_avg, bhp, /*iterate=*/ true, well_potentials, deferred_logger); for (double& value : well_potentials) { // make the value a little safer in case the BHP limits are default ones // TODO: a better way should be a better rescaling based on the investigation of the VFP table. @@ -2445,7 +2492,7 @@ namespace Opm } } - well_potentials = computeWellPotentialWithTHP(ebosSimulator, bhp, well_potentials, deferred_logger); + well_potentials = computeWellPotentialWithTHP(ebosSimulator, B_avg, bhp, well_potentials, deferred_logger); } } @@ -2861,7 +2908,7 @@ namespace Opm template void StandardWellV:: - wellTestingPhysical(Simulator& ebos_simulator, const std::vector& B_avg, + wellTestingPhysical(const Simulator& ebos_simulator, const std::vector& B_avg, const double /* simulation_time */ , const int /* report_step */, WellState& well_state, WellTestState& welltest_state, Opm::DeferredLogger& deferred_logger) { @@ -2884,7 +2931,7 @@ namespace Opm // we should be able to provide a better initialization calculateExplicitQuantities(ebos_simulator, well_state_copy, deferred_logger); - updateWellOperability(ebos_simulator, well_state_copy, deferred_logger); + updateWellOperability(ebos_simulator, B_avg, well_state_copy, deferred_logger); if ( !this->isOperable() ) { const std::string msg = " well " + name() + " is not operable during well testing for physical reason"; @@ -2892,7 +2939,7 @@ namespace Opm return; } - updateWellStateWithTarget(ebos_simulator, well_state_copy, deferred_logger); + updateWellStateWithTarget(ebos_simulator, B_avg, well_state_copy, deferred_logger); calculateExplicitQuantities(ebos_simulator, well_state_copy, deferred_logger); updatePrimaryVariables(well_state_copy, deferred_logger); diff --git a/opm/simulators/wells/StandardWell_impl.hpp b/opm/simulators/wells/StandardWell_impl.hpp index cd0ff3df7..eb38e4f76 100644 --- a/opm/simulators/wells/StandardWell_impl.hpp +++ b/opm/simulators/wells/StandardWell_impl.hpp @@ -444,13 +444,13 @@ namespace Opm void StandardWell:: assembleWellEq(const Simulator& ebosSimulator, - const std::vector& /* B_avg */, + const std::vector& B_avg, const double dt, WellState& well_state, Opm::DeferredLogger& deferred_logger) { - checkWellOperability(ebosSimulator, well_state, deferred_logger); + checkWellOperability(ebosSimulator, B_avg, well_state, deferred_logger); if (!this->isOperable()) return; @@ -1117,8 +1117,9 @@ namespace Opm void StandardWell:: updateWellStateWithTarget(const Simulator& ebos_simulator, + const std::vector& B_avg, WellState& well_state, - Opm::DeferredLogger& deferred_logger) const + Opm::DeferredLogger& deferred_logger) { // number of phases @@ -1144,7 +1145,7 @@ namespace Opm case THP: { // when a well can not work under THP target, it switches to BHP control if (this->operability_status_.isOperableUnderTHPLimit() ) { - updateWellStateWithTHPTargetIPR(ebos_simulator, well_state, deferred_logger); + updateWellStateWithTHPTargetIPR(ebos_simulator, B_avg, well_state, deferred_logger); } else { // go to BHP limit assert(this->operability_status_.isOperableUnderBHPLimit() ); @@ -1322,6 +1323,7 @@ namespace Opm void StandardWell:: checkWellOperability(const Simulator& ebos_simulator, + const std::vector& B_avg, const WellState& well_state, Opm::DeferredLogger& deferred_logger ) @@ -1337,7 +1339,7 @@ namespace Opm const bool old_well_operable = this->operability_status_.isOperable(); - updateWellOperability(ebos_simulator, well_state, deferred_logger); + updateWellOperability(ebos_simulator, B_avg, well_state, deferred_logger); const bool well_operable = this->operability_status_.isOperable(); @@ -1356,6 +1358,7 @@ namespace Opm void StandardWell:: updateWellOperability(const Simulator& ebos_simulator, + const std::vector& B_avg, const WellState& /* well_state */, Opm::DeferredLogger& deferred_logger ) @@ -1365,7 +1368,7 @@ namespace Opm updateIPR(ebos_simulator, deferred_logger); // checking the BHP limit related - checkOperabilityUnderBHPLimitProducer(ebos_simulator, deferred_logger); + checkOperabilityUnderBHPLimitProducer(ebos_simulator, B_avg, deferred_logger); // checking whether the well can operate under the THP constraints. if (this->wellHasTHPConstraints()) { @@ -1380,7 +1383,9 @@ namespace Opm template void StandardWell:: - checkOperabilityUnderBHPLimitProducer(const Simulator& ebos_simulator, Opm::DeferredLogger& deferred_logger) + checkOperabilityUnderBHPLimitProducer(const Simulator& ebos_simulator, + const std::vector& B_avg, + Opm::DeferredLogger& deferred_logger) { const double bhp_limit = mostStrictBhpFromBhpLimits(deferred_logger); // Crude but works: default is one atmosphere. @@ -1404,7 +1409,7 @@ namespace Opm // option 2: stick with the above IPR curve // we use IPR here std::vector well_rates_bhp_limit; - computeWellRatesWithBhp(ebos_simulator, bhp_limit, well_rates_bhp_limit, deferred_logger); + computeWellRatesWithBhp(ebos_simulator, B_avg, bhp_limit, /*iterate=*/ false, well_rates_bhp_limit, deferred_logger); const double thp = calculateThpFromBhp(well_rates_bhp_limit, bhp_limit, deferred_logger); const double thp_limit = this->getTHPConstraint(deferred_logger); @@ -1510,7 +1515,7 @@ namespace Opm { const double bhp = well_state.bhp()[index_of_well_]; std::vector well_rates; - computeWellRatesWithBhp(ebos_simulator, bhp, well_rates, deferred_logger); + computeWellRatesWithBhp(ebos_simulator, bhp, /*iterate=*/ false, well_rates, deferred_logger); const double sign = (well_type_ == PRODUCER) ? -1. : 1.; const double threshold = sign * std::numeric_limits::min(); @@ -1553,11 +1558,13 @@ namespace Opm void StandardWell:: updateWellStateWithTHPTargetIPR(const Simulator& ebos_simulator, + const std::vector& B_avg, WellState& well_state, - Opm::DeferredLogger& deferred_logger) const + Opm::DeferredLogger& deferred_logger) { if (well_type_ == PRODUCER) { updateWellStateWithTHPTargetIPRProducer(ebos_simulator, + B_avg, well_state, deferred_logger); } @@ -1577,8 +1584,9 @@ namespace Opm void StandardWell:: updateWellStateWithTHPTargetIPRProducer(const Simulator& ebos_simulator, + const std::vector& B_avg, WellState& well_state, - Opm::DeferredLogger& deferred_logger) const + Opm::DeferredLogger& deferred_logger) { well_state.thp()[index_of_well_] = this->getTHPConstraint(deferred_logger); @@ -1594,7 +1602,7 @@ namespace Opm initPrimaryVariablesEvaluation(); std::vector rates; - computeWellRatesWithBhp(ebos_simulator, bhp, rates, deferred_logger); + computeWellRatesWithBhp(ebos_simulator, B_avg, bhp, /*iterate=*/ false, rates, deferred_logger); // TODO: double checke the obtained rates // this is another places we might obtain negative rates @@ -1617,7 +1625,7 @@ namespace Opm calculateBHPWithTHPTargetIPR(Opm::DeferredLogger& deferred_logger) const { const double thp_target = this->getTHPConstraint(deferred_logger); - const double thp_control_index = this->getTHPControlIndex(); + const double thp_control_index = this->getControlIndex(THP); const int thp_table_id = well_controls_iget_vfp(well_controls_, thp_control_index); const double alq = well_controls_iget_alq(well_controls_, thp_control_index); @@ -2157,12 +2165,43 @@ namespace Opm void StandardWell:: computeWellRatesWithBhp(const Simulator& ebosSimulator, - const EvalWell& bhp, + const std::vector& B_avg, + const double& bhp, + const bool iterate, std::vector& well_flux, - Opm::DeferredLogger& deferred_logger) const + Opm::DeferredLogger& deferred_logger) { + + const int np = number_of_phases_; well_flux.resize(np, 0.0); + WellControls* wc = well_controls_; + const int bhp_index = Base::getControlIndex(BHP); + const double orig_bhp = well_controls_iget_target(wc, bhp_index); + const auto orig_current = well_controls_get_current(wc); + + well_controls_iset_target(wc, bhp_index, bhp); + well_controls_set_current(wc, bhp_index); + + // iterate to get a more accurate well density + if (iterate) { + // create a copy of the well_state to use. If the operability checking is sucessful, we use this one + // to replace the original one + WellState well_state_copy = ebosSimulator.problem().wellModel().wellState(); + well_state_copy.currentControls()[index_of_well_] = bhp_index; + + bool converged = this->solveWellEqUntilConverged(ebosSimulator, B_avg, well_state_copy, deferred_logger); + + if (!converged) { + const std::string msg = " well " + name() + " did not get converged during well potential calculations " + "returning zero values for the potential"; + deferred_logger.debug(msg); + return; + } + updatePrimaryVariables(well_state_copy, deferred_logger); + computeWellConnectionPressures(ebosSimulator, well_state_copy); + initPrimaryVariablesEvaluation(); + } const bool allow_cf = getAllowCrossFlow(); @@ -2185,6 +2224,11 @@ namespace Opm well_flux[ebosCompIdxToFlowCompIdx(p)] += cq_s[p].value(); } } + + + // reset bhp limit + well_controls_iset_target(wc, bhp_index, orig_bhp); + well_controls_set_current(wc, orig_current); } @@ -2195,9 +2239,10 @@ namespace Opm std::vector StandardWell:: computeWellPotentialWithTHP(const Simulator& ebosSimulator, + const std::vector& B_avg, const double initial_bhp, // bhp from BHP constraints const std::vector& initial_potential, - Opm::DeferredLogger& deferred_logger) const + Opm::DeferredLogger& deferred_logger) { // TODO: pay attention to the situation that finally the potential is calculated based on the bhp control // TODO: should we consider the bhp constraints during the iterative process? @@ -2261,7 +2306,7 @@ namespace Opm converged = std::abs(old_bhp - bhp) < bhp_tolerance; - computeWellRatesWithBhp(ebosSimulator, bhp, potentials, deferred_logger); + computeWellRatesWithBhp(ebosSimulator, B_avg, bhp, /*iterate=*/ true, potentials, deferred_logger); // checking whether the potentials have valid values for (const double value : potentials) { @@ -2303,6 +2348,7 @@ namespace Opm void StandardWell:: computeWellPotentials(const Simulator& ebosSimulator, + const std::vector& B_avg, const WellState& well_state, std::vector& well_potentials, Opm::DeferredLogger& deferred_logger) // const @@ -2315,6 +2361,7 @@ namespace Opm // TODO: for computeWellPotentials, no derivative is required actually initPrimaryVariablesEvaluation(); + const int np = number_of_phases_; well_potentials.resize(np, 0.0); @@ -2324,7 +2371,7 @@ namespace Opm // does the well have a THP related constraint? if ( !wellHasTHPConstraints() ) { assert(std::abs(bhp) != std::numeric_limits::max()); - computeWellRatesWithBhp(ebosSimulator, bhp, well_potentials, deferred_logger); + computeWellRatesWithBhp(ebosSimulator, B_avg, bhp, /*iterate=*/ true, well_potentials, deferred_logger); } else { // the well has a THP related constraint // checking whether a well is newly added, it only happens at the beginning of the report step @@ -2336,7 +2383,7 @@ namespace Opm } } else { // We need to generate a reasonable rates to start the iteration process - computeWellRatesWithBhp(ebosSimulator, bhp, well_potentials, deferred_logger); + computeWellRatesWithBhp(ebosSimulator, B_avg, bhp, /*iterate=*/ true, well_potentials, deferred_logger); for (double& value : well_potentials) { // make the value a little safer in case the BHP limits are default ones // TODO: a better way should be a better rescaling based on the investigation of the VFP table. @@ -2345,7 +2392,7 @@ namespace Opm } } - well_potentials = computeWellPotentialWithTHP(ebosSimulator, bhp, well_potentials, deferred_logger); + well_potentials = computeWellPotentialWithTHP(ebosSimulator, B_avg, bhp, well_potentials, deferred_logger); } } @@ -2754,7 +2801,7 @@ namespace Opm template void StandardWell:: - wellTestingPhysical(Simulator& ebos_simulator, const std::vector& B_avg, + wellTestingPhysical(const Simulator& ebos_simulator, const std::vector& B_avg, const double /* simulation_time */, const int /* report_step */, WellState& well_state, WellTestState& welltest_state, Opm::DeferredLogger& deferred_logger) @@ -2778,7 +2825,7 @@ namespace Opm // we should be able to provide a better initialization calculateExplicitQuantities(ebos_simulator, well_state_copy, deferred_logger); - updateWellOperability(ebos_simulator, well_state_copy, deferred_logger); + updateWellOperability(ebos_simulator, B_avg, well_state_copy, deferred_logger); if ( !this->isOperable() ) { const std::string msg = " well " + name() + " is not operable during well testing for physical reason"; @@ -2786,7 +2833,7 @@ namespace Opm return; } - updateWellStateWithTarget(ebos_simulator, well_state_copy, deferred_logger); + updateWellStateWithTarget(ebos_simulator, B_avg, well_state_copy, deferred_logger); calculateExplicitQuantities(ebos_simulator, well_state_copy, deferred_logger); updatePrimaryVariables(well_state_copy, deferred_logger); diff --git a/opm/simulators/wells/WellInterface.hpp b/opm/simulators/wells/WellInterface.hpp index 0fb501729..616c1f3cb 100644 --- a/opm/simulators/wells/WellInterface.hpp +++ b/opm/simulators/wells/WellInterface.hpp @@ -182,15 +182,18 @@ namespace Opm // TODO: before we decide to put more information under mutable, this function is not const virtual void computeWellPotentials(const Simulator& ebosSimulator, + const std::vector& B_avg, const WellState& well_state, std::vector& well_potentials, Opm::DeferredLogger& deferred_logger) = 0; virtual void updateWellStateWithTarget(const Simulator& ebos_simulator, + const std::vector& B_avg, WellState& well_state, - Opm::DeferredLogger& deferred_logger) const = 0; + Opm::DeferredLogger& deferred_logger) = 0; void updateWellControl(const Simulator& ebos_simulator, + const std::vector& B_avg, WellState& well_state, Opm::DeferredLogger& deferred_logger) /* const */; @@ -235,7 +238,7 @@ namespace Opm // TODO: theoretically, it should be a const function // Simulator is not const is because that assembleWellEq is non-const Simulator - void wellTesting(Simulator& simulator, const std::vector& Bavg, + void wellTesting(const Simulator& simulator, const std::vector& B_avg, const double simulation_time, const int report_step, const WellTestConfig::Reason testing_reason, /* const */ WellState& well_state, WellTestState& welltest_state, @@ -243,7 +246,10 @@ namespace Opm void updatePerforatedCell(std::vector& is_cell_perforated); - virtual void checkWellOperability(const Simulator& ebos_simulator, const WellState& well_state, Opm::DeferredLogger& deferred_logger) = 0; + virtual void checkWellOperability(const Simulator& ebos_simulator, + const std::vector& B_avg, + const WellState& well_state, + Opm::DeferredLogger& deferred_logger) = 0; // whether the well is operable bool isOperable() const; @@ -354,7 +360,7 @@ namespace Opm double getTHPConstraint(Opm::DeferredLogger& deferred_logger) const; - int getTHPControlIndex() const; + int getControlIndex(const WellControlType& type) const; // Component fractions for each phase for the well const std::vector& compFrac() const; @@ -385,11 +391,11 @@ namespace Opm OperabilityStatus operability_status_; - void wellTestingEconomic(Simulator& simulator, const std::vector& B_avg, + void wellTestingEconomic(const Simulator& simulator, const std::vector& B_avg, const double simulation_time, const int report_step, const WellState& well_state, WellTestState& welltest_state, Opm::DeferredLogger& deferred_logger); - virtual void wellTestingPhysical(Simulator& simulator, const std::vector& B_avg, + virtual void wellTestingPhysical(const Simulator& simulator, const std::vector& B_avg, const double simulation_time, const int report_step, WellState& well_state, WellTestState& welltest_state, Opm::DeferredLogger& deferred_logger) = 0; @@ -405,17 +411,25 @@ namespace Opm WellTestState& well_test_state, Opm::DeferredLogger& deferred_logger) const; - void solveWellForTesting(Simulator& ebosSimulator, WellState& well_state, + void solveWellForTesting(const Simulator& ebosSimulator, WellState& well_state, const std::vector& B_avg, Opm::DeferredLogger& deferred_logger); - bool solveWellEqUntilConverged(Simulator& ebosSimulator, + bool solveWellEqUntilConverged(const Simulator& ebosSimulator, const std::vector& B_avg, WellState& well_state, Opm::DeferredLogger& deferred_logger); void scaleProductivityIndex(const int perfIdx, double& productivity_index, const bool new_well, Opm::DeferredLogger& deferred_logger); + virtual void computeWellRatesWithBhp(const Simulator& ebosSimulator, + const std::vector& B_avg, + const double& bhp, + const bool iterate, + std::vector& well_flux, + Opm::DeferredLogger& deferred_logger) = 0; + + // count the number of times an output log message is created in the productivity // index calculations int well_productivity_index_logger_counter_; diff --git a/opm/simulators/wells/WellInterface_impl.hpp b/opm/simulators/wells/WellInterface_impl.hpp index 65b35360a..da6d9379f 100644 --- a/opm/simulators/wells/WellInterface_impl.hpp +++ b/opm/simulators/wells/WellInterface_impl.hpp @@ -380,7 +380,7 @@ namespace Opm WellInterface:: wellHasTHPConstraints() const { - return getTHPControlIndex() >= 0; + return getControlIndex(THP) >= 0; } @@ -391,7 +391,7 @@ namespace Opm WellInterface:: getTHPConstraint(Opm::DeferredLogger& deferred_logger) const { - const int thp_control_index = getTHPControlIndex(); + const int thp_control_index = getControlIndex(THP); if (thp_control_index < 0) { OPM_DEFLOG_THROW(std::runtime_error, " there is no THP constraint/limit for well " << name() @@ -407,11 +407,11 @@ namespace Opm template int WellInterface:: - getTHPControlIndex() const + getControlIndex(const WellControlType& type) const { const int nwc = well_controls_get_num(well_controls_); for (int ctrl_index = 0; ctrl_index < nwc; ++ctrl_index) { - if (well_controls_iget_type(well_controls_, ctrl_index) == THP) { + if (well_controls_iget_type(well_controls_, ctrl_index) == type) { return ctrl_index; } } @@ -426,6 +426,7 @@ namespace Opm void WellInterface:: updateWellControl(const Simulator& ebos_simulator, + const std::vector& B_avg, WellState& well_state, Opm::DeferredLogger& deferred_logger) /* const */ { @@ -504,7 +505,7 @@ namespace Opm } if (updated_control_index != old_control_index) { // || well_collection_->groupControlActive()) { - updateWellStateWithTarget(ebos_simulator, well_state, deferred_logger); + updateWellStateWithTarget(ebos_simulator, B_avg, well_state, deferred_logger); updatePrimaryVariables(well_state, deferred_logger); } } @@ -934,7 +935,7 @@ namespace Opm template void WellInterface:: - wellTesting(Simulator& simulator, const std::vector& B_avg, + wellTesting(const Simulator& simulator, const std::vector& B_avg, const double simulation_time, const int report_step, const WellTestConfig::Reason testing_reason, /* const */ WellState& well_state, @@ -959,7 +960,7 @@ namespace Opm template void WellInterface:: - wellTestingEconomic(Simulator& simulator, const std::vector& B_avg, + wellTestingEconomic(const Simulator& simulator, const std::vector& B_avg, const double simulation_time, const int report_step, const WellState& well_state, WellTestState& welltest_state, Opm::DeferredLogger& deferred_logger) { @@ -1166,7 +1167,7 @@ namespace Opm template bool WellInterface:: - solveWellEqUntilConverged(Simulator& ebosSimulator, + solveWellEqUntilConverged(const Simulator& ebosSimulator, const std::vector& B_avg, WellState& well_state, Opm::DeferredLogger& deferred_logger) @@ -1189,7 +1190,7 @@ namespace Opm ++it; solveEqAndUpdateWellState(well_state, deferred_logger); - updateWellControl(ebosSimulator, well_state, deferred_logger); + updateWellControl(ebosSimulator, B_avg, well_state, deferred_logger); initPrimaryVariablesEvaluation(); } while (it < max_iter); @@ -1238,7 +1239,7 @@ namespace Opm template void WellInterface:: - solveWellForTesting(Simulator& ebosSimulator, WellState& well_state, + solveWellForTesting(const Simulator& ebosSimulator, WellState& well_state, const std::vector& B_avg, Opm::DeferredLogger& deferred_logger) {