diff --git a/CMakeLists_files.cmake b/CMakeLists_files.cmake index a6c9947d9..8a34d6a26 100644 --- a/CMakeLists_files.cmake +++ b/CMakeLists_files.cmake @@ -70,6 +70,7 @@ list (APPEND TEST_SOURCE_FILES tests/test_stoppedwells.cpp tests/test_relpermdiagnostics.cpp tests/test_norne_pvt.cpp + tests/test_wellstatefullyimplicitblackoil.cpp ) if(MPI_FOUND) diff --git a/opm/autodiff/BlackoilWellModel.hpp b/opm/autodiff/BlackoilWellModel.hpp index 9fb5b71b1..57411760a 100644 --- a/opm/autodiff/BlackoilWellModel.hpp +++ b/opm/autodiff/BlackoilWellModel.hpp @@ -423,9 +423,15 @@ namespace Opm { // convert well data from opm-common to well state from opm-core void wellsToState( const data::Wells& wells, - PhaseUsage phases, - WellStateFullyImplicitBlackoil& state ); + const PhaseUsage& phases, + const bool handle_ms_well, + const int report_step, + WellStateFullyImplicitBlackoil& state ) const; + // whether there exists any multisegment well open on this process + bool anyMSWellOpenLocal(const Wells* wells, const int report_step) const; + + const Well* getWellEcl(const std::string& well_name) const; }; diff --git a/opm/autodiff/BlackoilWellModel_impl.hpp b/opm/autodiff/BlackoilWellModel_impl.hpp index 5caadc3e8..4fb5050ef 100644 --- a/opm/autodiff/BlackoilWellModel_impl.hpp +++ b/opm/autodiff/BlackoilWellModel_impl.hpp @@ -259,16 +259,8 @@ namespace Opm { well_state_.init(wells(), cellPressures, schedule(), wells_ecl_, timeStepIdx, &previous_well_state_, phase_usage_); // handling MS well related - if (param_.use_multisegment_well_) { // if we use MultisegmentWell model - for (const auto& well : wells_ecl_) { - // TODO: this is acutally not very accurate, because sometimes a deck just claims a MS well - // while keep the well shut. More accurately, we should check if the well exisits in the Wells - // structure here - if (well->isMultiSegment(timeStepIdx) ) { // there is one well is MS well - well_state_.initWellStateMSWell(wells(), wells_ecl_, timeStepIdx, phase_usage_, previous_well_state_); - break; - } - } + if (param_.use_multisegment_well_&& anyMSWellOpenLocal(wells(), timeStepIdx)) { // if we use MultisegmentWell model + well_state_.initWellStateMSWell(wells(), wells_ecl_, timeStepIdx, phase_usage_, &previous_well_state_); } // update the previous well state. This is used to restart failed steps. @@ -480,13 +472,16 @@ namespace Opm { initFromRestartFile(const RestartValue& restartValues) { const auto& defunctWellNames = ebosSimulator_.vanguard().defunctWellNames(); + + // The restart step value is used to identify wells present at the given + // time step. Wells that are added at the same time step as RESTART is initiated + // will not be present in a restart file. Use the previous time step to retrieve + // wells that have information written to the restart file. + const int report_step = std::max(eclState().getInitConfig().getRestartStep() - 1, 0); + WellsManager wellsmanager(eclState(), schedule(), - // The restart step value is used to identify wells present at the given - // time step. Wells that are added at the same time step as RESTART is initiated - // will not be present in a restart file. Use the previous time step to retrieve - // wells that have information written to the restart file. - std::max(eclState().getInitConfig().getRestartStep() - 1, 0), + report_step, Opm::UgGridHelpers::numCells(grid()), Opm::UgGridHelpers::globalCell(grid()), Opm::UgGridHelpers::cartDims(grid()), @@ -502,8 +497,9 @@ namespace Opm { if (nw > 0) { const auto phaseUsage = phaseUsageFromDeck(eclState()); const size_t numCells = Opm::UgGridHelpers::numCells(grid()); - well_state_.resize(wells, schedule(), numCells, phaseUsage); // Resize for restart step - wellsToState(restartValues.wells, phaseUsage, well_state_); + const bool handle_ms_well = (param_.use_multisegment_well_ && anyMSWellOpenLocal(wells, report_step)); + well_state_.resize(wells, wells_ecl_, schedule(), handle_ms_well, report_step, numCells, phaseUsage); // Resize for restart step + wellsToState(restartValues.wells, phaseUsage, handle_ms_well, report_step, well_state_); previous_well_state_ = well_state_; } initial_step_ = false; @@ -1667,8 +1663,11 @@ namespace Opm { void BlackoilWellModel:: wellsToState( const data::Wells& wells, - PhaseUsage phases, - WellStateFullyImplicitBlackoil& state ) { + const PhaseUsage& phases, + const bool handle_ms_well, + const int report_step, + WellStateFullyImplicitBlackoil& state ) const + { using rt = data::Rates::opt; const auto np = phases.num_phases; @@ -1724,7 +1723,80 @@ namespace Opm { } ++local_comp_index; } + + if (handle_ms_well && !well.segments.empty()) { + // we need the well_ecl_ information + const std::string& well_name = wm.first; + const Well* well_ecl = getWellEcl(well_name); + assert(well_ecl); + + const WellSegments& segment_set = well_ecl->getWellSegments(report_step); + + const int top_segment_index = state.topSegmentIndex(well_index); + const auto& segments = well.segments; + + // \Note: eventually we need to hanlde the situations that some segments are shut + assert(segment_set.size() == segments.size()); + + for (const auto& segment : segments) { + const int segment_index = segment_set.segmentNumberToIndex(segment.first); + + // recovering segment rates and pressure from the restart values + state.segPress()[top_segment_index + segment_index] = segment.second.pressure; + + const auto& segment_rates = segment.second.rates; + for (int p = 0; p < np; ++p) { + state.segRates()[(top_segment_index + segment_index) * np + p] = segment_rates.get(phs[p]); + } + } + } } } + + + + + template + bool + BlackoilWellModel:: + anyMSWellOpenLocal(const Wells* wells, const int report_step) const + { + bool any_ms_well_open = false; + + const int nw = wells->number_of_wells; + for (int w = 0; w < nw; ++w) { + const std::string well_name = std::string(wells->name[w]); + + const Well* well_ecl = getWellEcl(well_name); + + if (well_ecl->isMultiSegment(report_step) ) { + any_ms_well_open = true; + break; + } + } + return any_ms_well_open; + } + + + + + + template + const Well* + BlackoilWellModel:: + getWellEcl(const std::string& well_name) const + { + // finding the iterator of the well in wells_ecl + auto well_ecl = std::find_if(wells_ecl_.begin(), + wells_ecl_.end(), + [&well_name](const Well* elem)->bool { + return elem->name() == well_name; + }); + + assert(well_ecl != wells_ecl_.end()); + + return *well_ecl; + } + } // namespace Opm diff --git a/opm/autodiff/WellStateFullyImplicitBlackoil.hpp b/opm/autodiff/WellStateFullyImplicitBlackoil.hpp index c68cbbd12..e20325659 100644 --- a/opm/autodiff/WellStateFullyImplicitBlackoil.hpp +++ b/opm/autodiff/WellStateFullyImplicitBlackoil.hpp @@ -25,8 +25,12 @@ #include #include #include + +#include #include + #include + #include #include #include @@ -276,20 +280,29 @@ namespace Opm { // we need to create a trival segment related values to avoid there will be some // multi-segment wells added later. - top_segment_index_.reserve(nw); + nseg_ = nw; + top_segment_index_.resize(nw); + seg_number_.resize(nw); for (int w = 0; w < nw; ++w) { - top_segment_index_.push_back(w); + top_segment_index_[w] = w; + seg_number_[w] = 1; // Top segment is segment #1 } segpress_ = bhp(); segrates_ = wellRates(); } } - void resize(const Wells* wells, const Schedule& schedule, std::size_t numCells, const PhaseUsage& pu) + + void resize(const Wells* wells, const std::vector& wells_ecl, const Schedule& schedule, + const bool handle_ms_well, const int report_step, const size_t numCells, + const PhaseUsage& pu) { const std::vector tmp(numCells, 0.0); // <- UGLY HACK to pass the size - const std::vector wells_ecl; init(wells, tmp, schedule, wells_ecl, 0, nullptr, pu); + + if (handle_ms_well) { + initWellStateMSWell(wells, wells_ecl, report_step, pu, nullptr); + } } /// Allocate and initialize if wells is non-null. Also tries @@ -448,9 +461,13 @@ namespace Opm { // we need to create a trival segment related values to avoid there will be some // multi-segment wells added later. + nseg_ = nw; + seg_number_.clear(); top_segment_index_.reserve(nw); + seg_number_.reserve(nw); for (int w = 0; w < nw; ++w) { top_segment_index_.push_back(w); + seg_number_.push_back(1); // Top segment is segment #1 } segpress_ = bhp(); segrates_ = wellRates(); @@ -575,6 +592,13 @@ namespace Opm } } assert(local_comp_index == this->wells_->well_connpos[ w + 1 ] - this->wells_->well_connpos[ w ]); + + const auto nseg = this->numSegments(w); + for (auto seg_ix = 0*nseg; seg_ix < nseg; ++seg_ix) { + const auto seg_no = this->segmentNumber(w, seg_ix); + well.segments[seg_no] = + this->reportSegmentResults(pu, w, seg_ix, seg_no); + } } return res; @@ -582,9 +606,8 @@ namespace Opm /// init the MS well related. - template void initWellStateMSWell(const Wells* wells, const std::vector& wells_ecl, - const int time_step, const PhaseUsage& pu, const PrevWellState& prev_well_state) + const int time_step, const PhaseUsage& pu, const WellStateFullyImplicitBlackoil* prev_well_state) { // still using the order in wells const int nw = wells->number_of_wells; @@ -598,6 +621,7 @@ namespace Opm segpress_.reserve(nw); segrates_.clear(); segrates_.reserve(nw * numPhases()); + seg_number_.clear(); nseg_ = 0; // in the init function, the well rates and perforation rates have been initialized or copied from prevState @@ -621,6 +645,7 @@ namespace Opm top_segment_index_.push_back(nseg_); if ( !well_ecl->isMultiSegment(time_step) ) { // not multi-segment well nseg_ += 1; + seg_number_.push_back(1); // Assign single segment (top) as number 1. segpress_.push_back(bhp()[w]); const int np = numPhases(); for (int p = 0; p < np; ++p) { @@ -634,6 +659,9 @@ namespace Opm const int well_nseg = segment_set.size(); const int nperf = completion_set.size(); nseg_ += well_nseg; + for (auto segID = 0*well_nseg; segID < well_nseg; ++segID) { + this->seg_number_.push_back(segment_set[segID].segmentNumber()); + } // we need to know for each segment, how many perforation it has and how many segments using it as outlet_segment // that is why I think we should use a well model to initialize the WellState here std::vector> segment_perforations(well_nseg); @@ -709,13 +737,13 @@ namespace Opm assert(int(segpress_.size()) == nseg_); assert(int(segrates_.size()) == nseg_ * numPhases() ); - if (!prev_well_state.wellMap().empty()) { + if (prev_well_state && !prev_well_state->wellMap().empty()) { // copying MS well related - const auto& end = prev_well_state.wellMap().end(); + const auto& end = prev_well_state->wellMap().end(); const int np = numPhases(); for (int w = 0; w < nw; ++w) { const std::string name( wells->name[w] ); - const auto& it = prev_well_state.wellMap().find( name ); + const auto& it = prev_well_state->wellMap().find( name ); if (it != end) { // the well is found in the prev_well_state // TODO: the well with same name can change a lot, like they might not have same number of segments @@ -723,7 +751,7 @@ namespace Opm // for now, we just copy them. const int old_index_well = (*it).second[0]; const int new_index_well = w; - const int old_top_segment_index = prev_well_state.topSegmentIndex(old_index_well); + const int old_top_segment_index = prev_well_state->topSegmentIndex(old_index_well); const int new_top_segmnet_index = topSegmentIndex(new_index_well); int number_of_segment = 0; // if it is the last well in list @@ -734,11 +762,11 @@ namespace Opm } for (int i = 0; i < number_of_segment * np; ++i) { - segrates_[new_top_segmnet_index * np + i] = prev_well_state.segRates()[old_top_segment_index * np + i]; + segrates_[new_top_segmnet_index * np + i] = prev_well_state->segRates()[old_top_segment_index * np + i]; } for (int i = 0; i < number_of_segment; ++i) { - segpress_[new_top_segmnet_index + i] = prev_well_state.segPress()[old_top_segment_index + i]; + segpress_[new_top_segmnet_index + i] = prev_well_state->segPress()[old_top_segment_index + i]; } } } @@ -931,6 +959,66 @@ namespace Opm // Well potentials std::vector well_potentials_; + /// Map segment index to segment number, mostly for MS wells. + /// + /// Segment number (one-based) of j-th segment in i-th well is + /// \code + /// const auto top = topSegmentIndex(i); + /// const auto seg_No = seg_number_[top + j]; + /// \end + std::vector seg_number_; + + ::Opm::data::Segment + reportSegmentResults(const PhaseUsage& pu, + const int well_id, + const int seg_ix, + const int seg_no) const + { + auto seg_res = ::Opm::data::Segment{}; + + const auto seg_dof = + this->topSegmentIndex(well_id) + seg_ix; + + const auto* rate = + &this->segRates()[seg_dof * this->numPhases()]; + + seg_res.pressure = this->segPress()[seg_dof]; + + if (pu.phase_used[Water]) { + seg_res.rates.set(data::Rates::opt::wat, + rate[pu.phase_pos[Water]]); + } + + if (pu.phase_used[Oil]) { + seg_res.rates.set(data::Rates::opt::oil, + rate[pu.phase_pos[Oil]]); + } + + if (pu.phase_used[Gas]) { + seg_res.rates.set(data::Rates::opt::gas, + rate[pu.phase_pos[Gas]]); + } + + seg_res.segNumber = seg_no; + + return seg_res; + } + + int numSegments(const int well_id) const + { + const auto topseg = this->topSegmentIndex(well_id); + + return (well_id + 1 == this->numWells()) // Last well? + ? (this->numSegment() - topseg) + : (this->topSegmentIndex(well_id + 1) - topseg); + } + + int segmentNumber(const int well_id, const int seg_id) const + { + const auto top_offset = this->topSegmentIndex(well_id); + + return this->seg_number_[top_offset + seg_id]; + } }; } // namespace Opm diff --git a/tests/test_wellstatefullyimplicitblackoil.cpp b/tests/test_wellstatefullyimplicitblackoil.cpp new file mode 100644 index 000000000..3694f100e --- /dev/null +++ b/tests/test_wellstatefullyimplicitblackoil.cpp @@ -0,0 +1,307 @@ +/* + Copyright 2018 Equinor ASA. + + This file is part of the Open Porous Media Project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . +*/ + +#include + +#define BOOST_TEST_MODULE WellStateFIBOTest + +#include + +#include + +#include +#include +#include +#include + +#include +#include +#include +#include + +#include + +#include +#include + +struct Setup +{ + Setup(const std::string& filename) + : Setup(Opm::Parser{}.parseFile(filename)) + {} + + Setup(const Opm::Deck& deck) + : es (deck) + , pu (Opm::phaseUsageFromDeck(es)) + , grid (es.getInputGrid()) + , sched(deck, es) + {} + + Opm::EclipseState es; + Opm::PhaseUsage pu; + Opm::GridManager grid; + Opm::Schedule sched; +}; + +namespace { + Opm::WellStateFullyImplicitBlackoil + buildWellState(const Setup& setup, const std::size_t timeStep) + { + auto state = Opm::WellStateFullyImplicitBlackoil{}; + + const auto cpress = + std::vector(setup.grid.c_grid()->number_of_cells, + 100.0*Opm::unit::barsa); + + const Opm::WellsManager wmgr{ + setup.es, setup.sched, timeStep, *setup.grid.c_grid() + }; + + state.init(wmgr.c_wells(), cpress, setup.sched, + setup.sched.getWells(timeStep), + timeStep, nullptr, setup.pu); + + state.initWellStateMSWell(wmgr.c_wells(), + setup.sched.getWells(timeStep), + timeStep, setup.pu, nullptr); + + return state; + } + + + void setSegPress(const std::vector& wells, + const std::size_t tstep, + Opm::WellStateFullyImplicitBlackoil& wstate) + { + const auto nWell = wells.size(); + + auto& segPress = wstate.segPress(); + + for (auto wellID = 0*nWell; wellID < nWell; ++wellID) { + const auto* well = wells[wellID]; + const auto topSegIx = wstate.topSegmentIndex(wellID); + const auto pressTop = 100.0 * wellID; + + auto* press = &segPress[topSegIx]; + + press[0] = pressTop; + + if (! well->isMultiSegment(tstep)) { + continue; + } + + const auto& segSet = well->getWellSegments(tstep); + const auto nSeg = segSet.size(); + + for (auto segID = 0*nSeg + 1; segID < nSeg; ++segID) { + // One-based numbering scheme for segments. + const auto segNo = segSet[segID].segmentNumber(); + press[segNo - 1] = pressTop + 1.0*(segNo - 1); + } + } + } + + + void setSegRates(const std::vector& wells, + const std::size_t tstep, + const Opm::PhaseUsage& pu, + Opm::WellStateFullyImplicitBlackoil& wstate) + { + const auto wat = pu.phase_used[Opm::BlackoilPhases::Aqua]; + const auto iw = wat ? pu.phase_pos[Opm::BlackoilPhases::Aqua] : -1; + + const auto oil = pu.phase_used[Opm::BlackoilPhases::Liquid]; + const auto io = oil ? pu.phase_pos[Opm::BlackoilPhases::Liquid] : -1; + + const auto gas = pu.phase_used[Opm::BlackoilPhases::Vapour]; + const auto ig = gas ? pu.phase_pos[Opm::BlackoilPhases::Vapour] : -1; + + const auto np = wstate.numPhases(); + + const auto nWell = wells.size(); + + auto& segRates = wstate.segRates(); + + for (auto wellID = 0*nWell; wellID < nWell; ++wellID) { + const auto* well = wells[wellID]; + const auto topSegIx = wstate.topSegmentIndex(wellID); + const auto rateTop = 1000.0 * wellID; + + if (wat) { segRates[np*topSegIx + iw] = rateTop; } + if (oil) { segRates[np*topSegIx + io] = rateTop; } + if (gas) { segRates[np*topSegIx + ig] = rateTop; } + + if (! well->isMultiSegment(tstep)) { + continue; + } + + const auto& segSet = well->getWellSegments(tstep); + const auto nSeg = segSet.size(); + + for (auto segID = 0*nSeg + 1; segID < nSeg; ++segID) { + // One-based numbering scheme for segments. + const auto segNo = segSet[segID].segmentNumber(); + + auto* rates = &segRates[(topSegIx + segNo - 1) * np]; + + if (wat) { rates[iw] = rateTop + 100.0*(segNo - 1); } + if (oil) { rates[io] = rateTop + 200.0*(segNo - 1); } + if (gas) { rates[ig] = rateTop + 400.0*(segNo - 1); } + } + } + } +} // Anonymous + +BOOST_AUTO_TEST_SUITE(Segment) + +// --------------------------------------------------------------------- + +BOOST_AUTO_TEST_CASE(Linearisation) +{ + const Setup setup{ "msw.data" }; + const auto tstep = std::size_t{0}; + + const auto wstate = buildWellState(setup, tstep); + + BOOST_CHECK_EQUAL(wstate.numSegment(), 6 + 1); + + const auto& wells = setup.sched.getWells(tstep); + BOOST_CHECK_EQUAL(wells.size(), 2); + + const auto prod01_first = wells[0]->name() == "PROD01"; + + BOOST_CHECK_EQUAL(wstate.topSegmentIndex(0), 0); + BOOST_CHECK_EQUAL(wstate.topSegmentIndex(1), + prod01_first ? 6 : 1); +} + +// --------------------------------------------------------------------- + +BOOST_AUTO_TEST_CASE(Pressure) +{ + const Setup setup{ "msw.data" }; + const auto tstep = std::size_t{0}; + + auto wstate = buildWellState(setup, tstep); + + const auto& wells = setup.sched.getWells(tstep); + const auto prod01_first = wells[0]->name() == "PROD01"; + + setSegPress(wells, tstep, wstate); + + const auto rpt = wstate.report(setup.pu, setup.grid.c_grid()->global_cell); + + { + const auto& xw = rpt.at("INJE01"); + + BOOST_CHECK_EQUAL(xw.segments.size(), 1); // Top Segment + + const auto& xseg = xw.segments.at(1); + + BOOST_CHECK_EQUAL(xseg.segNumber, 1); + BOOST_CHECK_CLOSE(xseg.pressure, prod01_first ? 100.0 : 0.0, 1.0e-10); + } + + { + const auto expect_nSeg = 6; + const auto& xw = rpt.at("PROD01"); + + BOOST_CHECK_EQUAL(xw.segments.size(), expect_nSeg); + + const auto pressTop = prod01_first ? 0.0 : 100.0; + + for (auto segID = 0; segID < expect_nSeg; ++segID) { + const auto& xseg = xw.segments.at(segID + 1); + + BOOST_CHECK_EQUAL(xseg.segNumber, segID + 1); + BOOST_CHECK_CLOSE(xseg.pressure, pressTop + 1.0*segID, 1.0e-10); + } + } +} + +// --------------------------------------------------------------------- + +BOOST_AUTO_TEST_CASE(Rates) +{ + const Setup setup{ "msw.data" }; + const auto tstep = std::size_t{0}; + + auto wstate = buildWellState(setup, tstep); + + const auto& wells = setup.sched.getWells(tstep); + const auto prod01_first = wells[0]->name() == "PROD01"; + + const auto& pu = setup.pu; + + setSegRates(wells, tstep, pu, wstate); + + const auto rpt = wstate.report(pu, setup.grid.c_grid()->global_cell); + + const auto wat = pu.phase_used[Opm::BlackoilPhases::Aqua]; + const auto oil = pu.phase_used[Opm::BlackoilPhases::Liquid]; + const auto gas = pu.phase_used[Opm::BlackoilPhases::Vapour]; + + BOOST_CHECK(wat && oil && gas); + + { + const auto rateTop = prod01_first ? 1000.0 : 0.0; + + const auto& xw = rpt.at("INJE01"); + + BOOST_CHECK_EQUAL(xw.segments.size(), 1); // Top Segment + + const auto& xseg = xw.segments.at(1); + + BOOST_CHECK_EQUAL(xseg.segNumber, 1); + BOOST_CHECK_CLOSE(xseg.rates.get(Opm::data::Rates::opt::wat), + rateTop, 1.0e-10); + + BOOST_CHECK_CLOSE(xseg.rates.get(Opm::data::Rates::opt::oil), + rateTop, 1.0e-10); + + BOOST_CHECK_CLOSE(xseg.rates.get(Opm::data::Rates::opt::gas), + rateTop, 1.0e-10); + } + + { + const auto expect_nSeg = 6; + const auto& xw = rpt.at("PROD01"); + + BOOST_CHECK_EQUAL(xw.segments.size(), expect_nSeg); + + const auto rateTop = prod01_first ? 0.0 : 1000.0; + + for (auto segNum = 1; segNum <= expect_nSeg; ++segNum) { + const auto& xseg = xw.segments.at(segNum); + + BOOST_CHECK_EQUAL(xseg.segNumber, segNum); + + BOOST_CHECK_CLOSE(xseg.rates.get(Opm::data::Rates::opt::wat), + rateTop + 100.0*(segNum - 1), 1.0e-10); + + BOOST_CHECK_CLOSE(xseg.rates.get(Opm::data::Rates::opt::oil), + rateTop + 200.0*(segNum - 1), 1.0e-10); + + BOOST_CHECK_CLOSE(xseg.rates.get(Opm::data::Rates::opt::gas), + rateTop + 400.0*(segNum - 1), 1.0e-10); + } + } +} + +BOOST_AUTO_TEST_SUITE_END()