diff --git a/opm/autodiff/BlackoilWellModel_impl.hpp b/opm/autodiff/BlackoilWellModel_impl.hpp index f7c3b319e..78dbbefad 100644 --- a/opm/autodiff/BlackoilWellModel_impl.hpp +++ b/opm/autodiff/BlackoilWellModel_impl.hpp @@ -172,7 +172,7 @@ namespace Opm { const double p = fs.pressure(FluidSystem::oilPhaseIdx).value(); cellPressures[cellIdx] = p; } - well_state_.init(wells(), cellPressures, &previous_well_state_, phase_usage_); + well_state_.init(wells(), cellPressures, wells_ecl_, timeStepIdx, &previous_well_state_, phase_usage_); // handling MS well related if (param_.use_multisegment_well_) { // if we use MultisegmentWell model @@ -364,7 +364,7 @@ namespace Opm { if (nw > 0) { auto phaseUsage = phaseUsageFromDeck(eclState()); size_t numCells = Opm::UgGridHelpers::numCells(grid()); - well_state_.resize(wells, numCells, phaseUsage); //Resize for restart step + well_state_.resize(wells, numCells); //Resize for restart step wellsToState(restartValues.wells, phaseUsage, well_state_); previous_well_state_ = well_state_; } @@ -938,18 +938,19 @@ namespace Opm { WellControls* wc = well->wellControls(); const int control = well_controls_get_current(wc); well_state_.currentControls()[w] = control; - // TODO: for VFP control, the perf_densities are still zero here, investigate better - // way to handle it later. - well->updateWellStateWithTarget(well_state_); - // The wells are not considered to be newly added - // for next time step - if (well_state_.isNewWell(w) ) { - well_state_.setNewWell(w, false); + if (well_state_.effectiveEventsHappen(w) ) { + well->updateWellStateWithTarget(well_state_); + } + + // there is no new well control change input within a report step, + // so next time step, the well does not consider to have effective events anymore + if (well_state_.effectiveEventsHappen(w) ) { + well_state_.setEffeciveEventHappen(w, false); } } // end of for (int w = 0; w < nw; ++w) - + updatePrimaryVariables(); } @@ -1185,8 +1186,10 @@ namespace Opm { // it will not change the control mode, only update the targets wellCollection().updateWellTargets(well_state_.wellRates()); + // TODO: we should only do the well is involved in the update group targets for (auto& well : well_container_) { well->updateWellStateWithTarget(well_state_); + well->updatePrimaryVariables(well_state_); } } } diff --git a/opm/autodiff/MultisegmentWell_impl.hpp b/opm/autodiff/MultisegmentWell_impl.hpp index 74bf8af38..9966ae49d 100644 --- a/opm/autodiff/MultisegmentWell_impl.hpp +++ b/opm/autodiff/MultisegmentWell_impl.hpp @@ -373,8 +373,6 @@ namespace Opm break; } // end of switch - - updatePrimaryVariables(well_state); } diff --git a/opm/autodiff/SimulatorFullyImplicitBlackoilOutput.hpp b/opm/autodiff/SimulatorFullyImplicitBlackoilOutput.hpp index 7d50da859..c45a65c54 100644 --- a/opm/autodiff/SimulatorFullyImplicitBlackoilOutput.hpp +++ b/opm/autodiff/SimulatorFullyImplicitBlackoilOutput.hpp @@ -460,7 +460,7 @@ namespace Opm const Wells* wells = wellsmanager.c_wells(); size_t numCells = Opm::UgGridHelpers::numCells(grid); - wellstate.resize(wells, numCells, phaseUsage ); //Resize for restart step + wellstate.resize(wells, numCells); //Resize for restart step auto restart_values = eclIO_->loadRestart(solution_keys, extra_keys); solutionToSim( restart_values, phaseUsage, simulatorstate ); diff --git a/opm/autodiff/StandardWell_impl.hpp b/opm/autodiff/StandardWell_impl.hpp index a94c56612..11e32f56b 100644 --- a/opm/autodiff/StandardWell_impl.hpp +++ b/opm/autodiff/StandardWell_impl.hpp @@ -1214,8 +1214,6 @@ namespace Opm break; } // end of switch - - updatePrimaryVariables(well_state); } @@ -1908,13 +1906,14 @@ namespace Opm } 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 - if ( !well_state.isNewWell(index_of_well_) ) { + if ( !well_state.effectiveEventsHappen(index_of_well_) ) { for (int p = 0; p < np; ++p) { // This is dangerous for new added well // since we are not handling the initialization correctly for now well_potentials[p] = well_state.wellRates()[index_of_well_ * np + p]; } - } else { + } else + { // We need to generate a reasonable rates to start the iteration process computeWellRatesWithBhp(ebosSimulator, bhp, well_potentials); for (double& value : well_potentials) { diff --git a/opm/autodiff/WellInterface_impl.hpp b/opm/autodiff/WellInterface_impl.hpp index f46f60241..d48379f64 100644 --- a/opm/autodiff/WellInterface_impl.hpp +++ b/opm/autodiff/WellInterface_impl.hpp @@ -438,6 +438,7 @@ namespace Opm if (updated_control_index != old_control_index) { // || well_collection_->groupControlActive()) { updateWellStateWithTarget(well_state); + updatePrimaryVariables(well_state); } } diff --git a/opm/autodiff/WellStateFullyImplicitBlackoil.hpp b/opm/autodiff/WellStateFullyImplicitBlackoil.hpp index 8f27b5d00..bdcd4c4df 100644 --- a/opm/autodiff/WellStateFullyImplicitBlackoil.hpp +++ b/opm/autodiff/WellStateFullyImplicitBlackoil.hpp @@ -62,7 +62,9 @@ namespace Opm /// Allocate and initialize if wells is non-null. Also tries /// to give useful initial values to the bhp(), wellRates() /// and perfPhaseRates() fields, depending on controls - void init(const Wells* wells, const std::vector& cellPressures, const WellStateFullyImplicitBlackoil* prevState, const PhaseUsage& pu) + void init(const Wells* wells, const std::vector& cellPressures, + const std::vector& wells_ecl, const int report_step, + const WellStateFullyImplicitBlackoil* prevState, const PhaseUsage& pu) { // call init on base class BaseType :: init(wells, cellPressures); @@ -84,15 +86,27 @@ namespace Opm well_dissolved_gas_rates_.resize(nw, 0.0); well_vaporized_oil_rates_.resize(nw, 0.0); - is_new_well_.resize(nw, true); - if (prevState && !prevState->wellMap().empty()) { - const auto& end = prevState->wellMap().end(); - for (int w = 0; w < nw; ++w) { - const auto& it = prevState->wellMap().find( wells->name[w]); - if (it != end) { - is_new_well_[w] = false; + // checking whether some effective well control happens + effective_events_happen_.resize(nw); + for (int w = 0; w name[w]); + for (; index_well_ecl < nw_wells_ecl; ++index_well_ecl) { + if (well_name == wells_ecl[index_well_ecl]->name()) { + break; } } + + // It should be able to find in wells_ecl. + if (index_well_ecl == nw_wells_ecl) { + OPM_THROW(std::logic_error, "Could not find well " << well_name << " in wells_ecl "); + } + + const Well* well_ecl = wells_ecl[index_well_ecl]; + // PRODUCTION_UPDATE, INJECTION_UPDATE, WELL_STATUS_CHANGE + // 16 + 32 + 128 + effective_events_happen_[w] = (well_ecl->hasEvent(176, report_step) ); } // Ensure that we start out with zero rates by default. @@ -117,9 +131,11 @@ namespace Opm } } - // Initialize current_controls_. - // The controls set in the Wells object are treated as defaults, - // and also used for initial values. + // set up the current well controls, if there are some effective events happened to the well, + // we use the control specified in the DECK + // if not, we will just use the control from the prevState + // the second part will be done when copying the data from prevState + // here, we initilaize all the current_controls_ based on the control mode from the DECK current_controls_.resize(nw); for (int w = 0; w < nw; ++w) { current_controls_[w] = well_controls_get_current(wells->ctrls[w]); @@ -136,7 +152,7 @@ namespace Opm typedef typename WellMapType :: const_iterator const_iterator; const_iterator end = prevState->wellMap().end(); for (int w = 0; w < nw; ++w) { - std::string name( wells->name[ w ] ); + const std::string name( wells->name[ w ] ); const_iterator it = prevState->wellMap().find( name ); if( it != end ) { @@ -149,6 +165,14 @@ namespace Opm // thp thp()[ newIndex ] = prevState->thp()[ oldIndex ]; + // if there is no effective control event happens to the well, we use the current_controls_ from prevState + if (!effective_events_happen_[w]) { + current_controls_[ newIndex ] = prevState->currentControls()[ oldIndex ]; + // also change the one in the WellControls + // TODO: checking if this is appropriate + well_controls_set_current(wells->ctrls[w], current_controls_[ newIndex ]); + } + // wellrates for( int i=0, idx=newIndex*np, oldidx=oldIndex*np; i tmp(numCells, 0.0); // <- UGLY HACK to pass the size - init(wells, tmp, nullptr, pu); + BaseType :: init(wells, tmp); + + const int np = wells->number_of_phases; + const int nw = wells->number_of_wells; + const int nperf = wells->well_connpos[nw]; + + well_reservoir_rates_.resize(nw * np, 0.0); + well_dissolved_gas_rates_.resize(nw, 0.0); + well_vaporized_oil_rates_.resize(nw, 0.0); + effective_events_happen_.resize(nw); + perfphaserates_.resize(nperf * np, 0.0); + current_controls_.resize(nw); + perfRateSolvent_.resize(nperf, 0.0); + top_segment_index_.resize(nw); + // TODO: segrates_ and segpress_ are not taken care of here } /// Allocate and initialize if wells is non-null. Also tries @@ -296,8 +334,6 @@ namespace Opm current_controls_[w] = well_controls_get_current(wells->ctrls[w]); } - is_new_well_.resize(nw, true); - perfRateSolvent_.clear(); perfRateSolvent_.resize(nperf, 0.0); @@ -312,9 +348,6 @@ namespace Opm const_iterator it = prevState.wellMap().find( name ); if( it != end ) { - // this is not a new added well - is_new_well_[w] = false; - const int oldIndex = (*it).second[ 0 ]; const int newIndex = w; @@ -707,13 +740,13 @@ namespace Opm } - bool isNewWell(const int w) const { - return is_new_well_[w]; + bool effectiveEventsHappen(const int w) const { + return effective_events_happen_[w]; } - void setNewWell(const int w, const bool is_new_well) { - is_new_well_[w] = is_new_well; + void setEffeciveEventHappen(const int w, const bool effective_events_happen) { + effective_events_happen_[w] = effective_events_happen; } @@ -802,11 +835,10 @@ namespace Opm // should be zero for injection wells std::vector well_vaporized_oil_rates_; - // marking whether the well is just added - // for newly added well, the current initialized rates from WellState - // will have very wrong compositions for production wells, will mostly cause - // problem with VFP interpolation - std::vector is_new_well_; + // some events happens to the well, like this well is a new well + // or new well control keywords happens + // \Note: for now, only WCON* keywords, and well status change is considered + std::vector effective_events_happen_; // MS well related // for StandardWell, the number of segments will be one