Merge pull request #1592 from bska/report-segment-res

WellState FIBO: Return Segment Vectors from report()
This commit is contained in:
Atgeirr Flø Rasmussen 2019-04-05 10:49:14 +02:00 committed by GitHub
commit a8d8fb86c3
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
5 changed files with 507 additions and 33 deletions

View File

@ -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)

View File

@ -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;
};

View File

@ -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<TypeTag>::
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<typename TypeTag>
bool
BlackoilWellModel<TypeTag>::
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<typename TypeTag>
const Well*
BlackoilWellModel<TypeTag>::
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

View File

@ -25,8 +25,12 @@
#include <opm/core/well_controls.h>
#include <opm/core/simulator/WellState.hpp>
#include <opm/core/props/BlackoilPhases.hpp>
#include <opm/parser/eclipse/EclipseState/Schedule/Schedule.hpp>
#include <opm/parser/eclipse/EclipseState/Schedule/Well/Well.hpp>
#include <opm/common/ErrorMacros.hpp>
#include <vector>
#include <cassert>
#include <string>
@ -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<const Well*>& wells_ecl, const Schedule& schedule,
const bool handle_ms_well, const int report_step, const size_t numCells,
const PhaseUsage& pu)
{
const std::vector<double> tmp(numCells, 0.0); // <- UGLY HACK to pass the size
const std::vector<const Well*> 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 <typename PrevWellState>
void initWellStateMSWell(const Wells* wells, const std::vector<const Well*>& 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<std::vector<int>> 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<double> 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<int> 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

View File

@ -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 <http://www.gnu.org/licenses/>.
*/
#include <config.h>
#define BOOST_TEST_MODULE WellStateFIBOTest
#include <opm/autodiff/WellStateFullyImplicitBlackoil.hpp>
#include <boost/test/unit_test.hpp>
#include <opm/parser/eclipse/Parser/Parser.hpp>
#include <opm/parser/eclipse/Parser/ParseContext.hpp>
#include <opm/parser/eclipse/EclipseState/Schedule/Schedule.hpp>
#include <opm/parser/eclipse/Units/Units.hpp>
#include <opm/core/props/BlackoilPhases.hpp>
#include <opm/core/props/phaseUsageFromDeck.hpp>
#include <opm/core/wells/WellsManager.hpp>
#include <opm/core/wells.h>
#include <opm/grid/GridManager.hpp>
#include <cstddef>
#include <string>
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<double>(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<const Opm::Well*>& 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<const Opm::Well*>& 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()