Merge pull request #3784 from bska/fix-ilbr4-ilbs-take-ii

Revise Branch Discovery Order for MSW Output
This commit is contained in:
Bård Skaflestad 2023-11-24 16:49:02 +01:00 committed by GitHub
commit 9a0ccfd853
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
2 changed files with 432 additions and 41 deletions

View File

@ -925,9 +925,9 @@ namespace {
/// 1, 2, 3, 4, 5, 6 -- Branch (1)
/// 11, 12, 13, 14, 15, 16 -- Branch (2)
/// 7, 8, 9, 10 -- Branch (3)
/// 17, 18, 19 -- Branch (4)
/// 20, 22, 23, 24 -- Branch (5)
/// 21, -- Branch (6)
/// 17, 18, 19 -- Branch (4)
///
/// +------------------------------------------------------------+
/// | |
@ -1007,6 +1007,49 @@ namespace {
void traverseStructure();
private:
/// Representation of a branch kick-off point
struct KickOffPoint
{
/// Branch start segment
int segment;
/// Branch outlet segment. Segment from which the branch
/// kicks off.
int outlet;
/// ID of branch starting at this kick-off point.
int branch;
};
/// Priority queue ordering operation
struct KickOffPointPriorityOrder
{
/// Priority comparison operation
///
/// \param[in] p1 Left-hand side of comparison
/// \param[in] p2 Right-hand side of comparison
///
/// \return Whether \p p1 has a lower priority than \p p2.
bool operator()(const KickOffPoint& p1,
const KickOffPoint& p2) const
{
// We use '>' here to have low outlet segment IDs move
// to the front of the priority queue. If multiple
// branches kick off from the same outlet segment, we
// sort these by their branch IDs whence lower-ID
// branches are created before higher-ID branches.
return std::tie(p1.outlet, p1.branch)
> std::tie(p2.outlet, p2.branch);
}
};
/// Priority queue container type for branch kick-off points.
using KickOffPointsQueue = std::priority_queue<
KickOffPoint,
std::vector<KickOffPoint>,
KickOffPointPriorityOrder
>;
/// Name of well being explored.
std::string_view well_{};
@ -1019,12 +1062,18 @@ namespace {
/// Callback routine for visiting a new branch.
NewBranchCallback runNewBranchCallback_{};
/// Segments from which to kick off searching new branches.
std::queue<int> kickOffSegments_{};
/// Points from which to kick off new branches.
KickOffPointsQueue kickOffPoints_{};
/// One-based segment number of currently visited segment.
int currentSegment_{};
/// Create a new branch from "top" kick-off point.
///
/// Invokes new branch callback and makes the "current" segment
/// be the kick-off segment.
void createNewBranch();
/// Find all segments on current branch, in order from heel to
/// toe.
///
@ -1042,20 +1091,6 @@ namespace {
void discoverNewBranches(const int outletSegment,
const std::vector<int>& children);
/// Enqueue new branch.
///
/// Will be visited later. Invokes new branch callback.
///
/// \param[in] branchId Branch number for new branch.
///
/// \param[in] kickOffSegment First segment on new branch.
///
/// \param[in] outletSegment Segment on branch from which the
/// new branch kicks off.
void discoverNewBranch(const int branchId,
const int kickOffSegment,
const int outletSegment);
/// Split child segments of current segment into groups
/// based on their associate branch number.
///
@ -1083,16 +1118,33 @@ namespace {
void Topology::traverseStructure()
{
this->kickOffSegments_.push(1);
while (! this->kickOffSegments_.empty()) {
this->currentSegment_ = this->kickOffSegments_.front();
this->kickOffSegments_.pop();
// Search starts at top segment: branch 1, segment 1.
this->kickOffPoints_.push(KickOffPoint { 1, 0, 1 });
while (! this->kickOffPoints_.empty()) {
this->createNewBranch();
this->buildCurrentBranch();
}
}
void Topology::createNewBranch()
{
const auto kickOff = this->kickOffPoints_.top();
this->kickOffPoints_.pop();
if (kickOff.branch != 1) {
// All branches other than the main branch (branch ID 1)
// create new branch objects. Run new branch callback for
// these.
this->runNewBranchCallback_(this->well_,
kickOff.branch,
kickOff.segment,
kickOff.outlet);
}
this->currentSegment_ = kickOff.segment;
}
void Topology::buildCurrentBranch()
{
while (true) {
@ -1121,23 +1173,14 @@ namespace {
const std::vector<int>& children)
{
for (const auto& child : children) {
const auto branch = this->segment(child).branchNumber();
this->discoverNewBranch(branch, child, outletSegment);
this->kickOffPoints_
.push(KickOffPoint {
child, outletSegment,
this->segment(child).branchNumber()
});
}
}
void Topology::discoverNewBranch(const int branchId,
const int kickOffSegment,
const int outletSegment)
{
this->runNewBranchCallback_(this->well_,
branchId,
kickOffSegment,
outletSegment);
this->kickOffSegments_.push(kickOffSegment);
}
std::pair<std::vector<int>, std::optional<int>>
Topology::characteriseChildSegments() const
{

View File

@ -243,6 +243,113 @@ END
)");
}
//------------------------------------------------------------------+
// Multi-lateral well with several ICDs and/or valves in a segment |
// structure of the form |
//------------------------------------------------------------------+
// |
// 8 9 10 11 12 |
// o-----o-----o------o------o------o (2) |
// 7 / \ |
// / \ |
// / 19 \ 20 26 |
// 1 2 3 / 4 5 6 o----o----o (4) |
// ---o---o---o-----o---o---o (1) |
// \ |
// 13 \ (7) (6) |
// \ o o |
// o 25 / 24 / |
// \ / / |
// 14 \ 15 / 16 17 / 18 |
// o-----o--------o-------o--------o (3) |
// \ |
// 21 \ 22 23 |
// o-------o--------o (5) |
// \ \ \ |
// 27 \ 28 \ 29 \ |
// o o o |
// (8) (9) (10) |
// |
//------------------------------------------------------------------+
// Branch ( 1): 1, 2, 3, 4, 5, 6 |
// Branch ( 2): 7, 8, 9, 10, 11, 12 |
// Branch ( 3): 13, 14, 15, 16, 17, 18 |
// Branch ( 4): 19, 20, 26 |
// Branch ( 5): 21, 22, 23 |
// Branch ( 6): 24 |
// Branch ( 7): 25 |
// Branch ( 8): 27 |
// Branch ( 9): 24 |
// Branch (10): 29 |
//------------------------------------------------------------------+
Opm::Deck multilaterals_with_icd_valve()
{
return Opm::Parser{}.parseString(R"(RUNSPEC
START
23 'NOV' 2023 /
DIMENS
10 10 3 /
OIL
GAS
WATER
DISGAS
VAPOIL
GRID
DXV
10*100.0 /
DYV
10*100.0 /
DZV
3*5.0 /
PERMX
300*100.0 /
COPY
PERMX PERMY /
PERMX PERMZ /
/
MULTIPLY
PERMZ 0.1 /
/
PORO
300*0.3 /
DEPTHZ
121*2000.0 /
SCHEDULE
WELSPECS
'MLP' 'G' 10 10 2002.5 'OIL' /
/
COMPDAT
'MLP' 10 10 3 3 'OPEN' 1* 123.4 /
/
WELSEGS
'MLP' 2002.5 0.0 1* 'INC' 'H--' /
--
2 6 1 1 0.1 0.1 0.2 0.01 /
7 12 2 3 0.1 0.1 0.2 0.01 /
13 18 3 2 0.1 0.1 0.2 0.01 /
19 20 4 10 0.1 0.1 0.2 0.01 /
26 26 4 20 0.1 0.1 0.2 0.01 /
21 23 5 15 0.1 0.1 0.2 0.01 /
24 24 6 17 0.1 0.1 0.2 0.01 /
25 25 7 15 0.1 0.1 0.2 0.01 /
27 27 8 21 0.1 0.1 0.2 0.01 /
28 28 9 22 0.1 0.1 0.2 0.01 /
29 29 10 23 0.1 0.1 0.2 0.01 /
/
COMPSEGS
'MLP' /
--
10 10 3 5 0.0 1.0 'Z' /
/
WCONPROD
'MLP' 'OPEN' 'ORAT' 321.0 4* 10.0 /
/
TSTEP
5*30 /
END
)");
}
} // Anonymous namespace
struct SimulationCase
@ -503,9 +610,9 @@ BOOST_AUTO_TEST_CASE (Declared_MSW_Data)
// 1, 2, 3, 4, 5, 6 -- Branch (1)
// 11, 12, 13, 14, 15, 16 -- Branch (2)
// 7, 8, 9, 10 -- Branch (3)
// 17, 18, 19 -- Branch (4)
// 20, 22, 23, 24 -- Branch (5)
// 21, -- Branch (6)
// 17, 18, 19 -- Branch (4)
//
BOOST_AUTO_TEST_CASE(Multilateral_Branches)
{
@ -540,7 +647,7 @@ BOOST_AUTO_TEST_CASE(Multilateral_Branches)
BOOST_CHECK_EQUAL(ilbs.size(), std::vector<int>::size_type{6});
const auto expect = std::vector {
11, 7, 20, 21, 17, 0,
11, 7, 17, 20, 21, 0,
};
BOOST_CHECK_EQUAL_COLLECTIONS(ilbs .begin(), ilbs .end(),
@ -593,7 +700,7 @@ BOOST_AUTO_TEST_CASE(Multilateral_Branches)
BOOST_CHECK_EQUAL(ilbr[VI::ILbr::NumBranchSegments], 3);
BOOST_CHECK_EQUAL(ilbr[VI::ILbr::FirstSegment], 17);
BOOST_CHECK_EQUAL(ilbr[VI::ILbr::LastSegment], 19);
BOOST_CHECK_EQUAL(ilbr[VI::ILbr::KickOffDiscoveryOffset], 5);
BOOST_CHECK_EQUAL(ilbr[VI::ILbr::KickOffDiscoveryOffset], 3);
}
// ILBR, branch 5
@ -604,7 +711,7 @@ BOOST_AUTO_TEST_CASE(Multilateral_Branches)
BOOST_CHECK_EQUAL(ilbr[VI::ILbr::NumBranchSegments], 4);
BOOST_CHECK_EQUAL(ilbr[VI::ILbr::FirstSegment], 20);
BOOST_CHECK_EQUAL(ilbr[VI::ILbr::LastSegment], 24);
BOOST_CHECK_EQUAL(ilbr[VI::ILbr::KickOffDiscoveryOffset], 3);
BOOST_CHECK_EQUAL(ilbr[VI::ILbr::KickOffDiscoveryOffset], 4);
}
// ILBR, branch 6
@ -615,7 +722,7 @@ BOOST_AUTO_TEST_CASE(Multilateral_Branches)
BOOST_CHECK_EQUAL(ilbr[VI::ILbr::NumBranchSegments], 1);
BOOST_CHECK_EQUAL(ilbr[VI::ILbr::FirstSegment], 21);
BOOST_CHECK_EQUAL(ilbr[VI::ILbr::LastSegment], 21);
BOOST_CHECK_EQUAL(ilbr[VI::ILbr::KickOffDiscoveryOffset], 4);
BOOST_CHECK_EQUAL(ilbr[VI::ILbr::KickOffDiscoveryOffset], 5);
}
}
@ -680,6 +787,247 @@ BOOST_AUTO_TEST_CASE(Multilateral_Segments_ISEG_0)
}
}
// The segments and branches must appear in the following order in the
// ILBS/ILBR output arrays.
//
// 1, 2, 3, 4, 5, 6 -- Branch ( 1)
// 13, 14, 15, 16, 17, 18 -- Branch ( 3)
// 7, 8, 9, 10, 11, 12 -- Branch ( 2)
// 19, 20, 26 -- Branch ( 4)
// 21, 22, 23 -- Branch ( 5)
// 25, -- Branch ( 7)
// 24, -- Branch ( 6)
// 27, -- Branch ( 8)
// 28, -- Branch ( 9)
// 29, -- Branch (10)
//
BOOST_AUTO_TEST_CASE(Multilateral_Branches_ICD_Valve)
{
const auto cse = SimulationCase { multilaterals_with_icd_valve() };
const auto& es = cse.es;
const auto& grid = cse.grid;
const auto& sched = cse.sched;
// Report Step 1: 2023-11-23 --> 2023-12-23
const auto rptStep = std::size_t {1};
const auto secs_elapsed = 30 * 86'400.0;
const auto ih = Opm::RestartIO::Helpers::
createInteHead(es, grid, sched, secs_elapsed,
rptStep, rptStep + 1, rptStep);
// No dynamic data needed for this test. We're checking the structure
// only.
const auto smry = Opm::SummaryState { Opm::TimeService::now() };
const auto xw = Opm::data::Wells {};
auto amswd = Opm::RestartIO::Helpers::AggregateMSWData {ih};
amswd.captureDeclaredMSWData(sched, rptStep, es.getUnits(),
ih, grid, smry, xw);
// ILBS--First segment on each branch other than branch 1. Ordered by
// discovery.
{
const auto& ilbs = amswd.getILBs();
// No WSEGDIMS => size = maximum branch number
BOOST_CHECK_EQUAL(ilbs.size(), std::vector<int>::size_type{10});
const auto expect = std::vector {
13, 7, 19, 21, 25, 24, 27, 28, 29, 0,
};
BOOST_CHECK_EQUAL_COLLECTIONS(ilbs .begin(), ilbs .end(),
expect.begin(), expect.end());
}
auto ilbrOffset = [&ih](const int branch)
{
return ih[VI::intehead::NILBRZ] * (branch - 1);
};
// ILBR, branch 1
{
const auto* ilbr = &amswd.getILBr()[ilbrOffset(1)];
BOOST_CHECK_EQUAL(ilbr[VI::ILbr::OutletSegment], 0);
BOOST_CHECK_EQUAL(ilbr[VI::ILbr::NumBranchSegments], 6);
BOOST_CHECK_EQUAL(ilbr[VI::ILbr::FirstSegment], 1);
BOOST_CHECK_EQUAL(ilbr[VI::ILbr::LastSegment], 6);
BOOST_CHECK_EQUAL(ilbr[VI::ILbr::KickOffDiscoveryOffset], 0);
}
// ILBR, branch 2
{
const auto* ilbr = &amswd.getILBr()[ilbrOffset(2)];
BOOST_CHECK_EQUAL(ilbr[VI::ILbr::OutletSegment], 3);
BOOST_CHECK_EQUAL(ilbr[VI::ILbr::NumBranchSegments], 6);
BOOST_CHECK_EQUAL(ilbr[VI::ILbr::FirstSegment], 7);
BOOST_CHECK_EQUAL(ilbr[VI::ILbr::LastSegment], 12);
BOOST_CHECK_EQUAL(ilbr[VI::ILbr::KickOffDiscoveryOffset], 2);
}
// ILBR, branch 3
{
const auto* ilbr = &amswd.getILBr()[ilbrOffset(3)];
BOOST_CHECK_EQUAL(ilbr[VI::ILbr::OutletSegment], 2);
BOOST_CHECK_EQUAL(ilbr[VI::ILbr::NumBranchSegments], 6);
BOOST_CHECK_EQUAL(ilbr[VI::ILbr::FirstSegment], 13);
BOOST_CHECK_EQUAL(ilbr[VI::ILbr::LastSegment], 18);
BOOST_CHECK_EQUAL(ilbr[VI::ILbr::KickOffDiscoveryOffset], 1);
}
// ILBR, branch 4
{
const auto* ilbr = &amswd.getILBr()[ilbrOffset(4)];
BOOST_CHECK_EQUAL(ilbr[VI::ILbr::OutletSegment], 10);
BOOST_CHECK_EQUAL(ilbr[VI::ILbr::NumBranchSegments], 3);
BOOST_CHECK_EQUAL(ilbr[VI::ILbr::FirstSegment], 19);
BOOST_CHECK_EQUAL(ilbr[VI::ILbr::LastSegment], 26);
BOOST_CHECK_EQUAL(ilbr[VI::ILbr::KickOffDiscoveryOffset], 3);
}
// ILBR, branch 5
{
const auto* ilbr = &amswd.getILBr()[ilbrOffset(5)];
BOOST_CHECK_EQUAL(ilbr[VI::ILbr::OutletSegment], 15);
BOOST_CHECK_EQUAL(ilbr[VI::ILbr::NumBranchSegments], 3);
BOOST_CHECK_EQUAL(ilbr[VI::ILbr::FirstSegment], 21);
BOOST_CHECK_EQUAL(ilbr[VI::ILbr::LastSegment], 23);
BOOST_CHECK_EQUAL(ilbr[VI::ILbr::KickOffDiscoveryOffset], 4);
}
// ILBR, branch 6
{
const auto* ilbr = &amswd.getILBr()[ilbrOffset(6)];
BOOST_CHECK_EQUAL(ilbr[VI::ILbr::OutletSegment], 17);
BOOST_CHECK_EQUAL(ilbr[VI::ILbr::NumBranchSegments], 1);
BOOST_CHECK_EQUAL(ilbr[VI::ILbr::FirstSegment], 24);
BOOST_CHECK_EQUAL(ilbr[VI::ILbr::LastSegment], 24);
BOOST_CHECK_EQUAL(ilbr[VI::ILbr::KickOffDiscoveryOffset], 6);
}
// ILBR, branch 7
{
const auto* ilbr = &amswd.getILBr()[ilbrOffset(7)];
BOOST_CHECK_EQUAL(ilbr[VI::ILbr::OutletSegment], 15);
BOOST_CHECK_EQUAL(ilbr[VI::ILbr::NumBranchSegments], 1);
BOOST_CHECK_EQUAL(ilbr[VI::ILbr::FirstSegment], 25);
BOOST_CHECK_EQUAL(ilbr[VI::ILbr::LastSegment], 25);
BOOST_CHECK_EQUAL(ilbr[VI::ILbr::KickOffDiscoveryOffset], 5);
}
// ILBR, branch 8
{
const auto* ilbr = &amswd.getILBr()[ilbrOffset(8)];
BOOST_CHECK_EQUAL(ilbr[VI::ILbr::OutletSegment], 21);
BOOST_CHECK_EQUAL(ilbr[VI::ILbr::NumBranchSegments], 1);
BOOST_CHECK_EQUAL(ilbr[VI::ILbr::FirstSegment], 27);
BOOST_CHECK_EQUAL(ilbr[VI::ILbr::LastSegment], 27);
BOOST_CHECK_EQUAL(ilbr[VI::ILbr::KickOffDiscoveryOffset], 7);
}
// ILBR, branch 9
{
const auto* ilbr = &amswd.getILBr()[ilbrOffset(9)];
BOOST_CHECK_EQUAL(ilbr[VI::ILbr::OutletSegment], 22);
BOOST_CHECK_EQUAL(ilbr[VI::ILbr::NumBranchSegments], 1);
BOOST_CHECK_EQUAL(ilbr[VI::ILbr::FirstSegment], 28);
BOOST_CHECK_EQUAL(ilbr[VI::ILbr::LastSegment], 28);
BOOST_CHECK_EQUAL(ilbr[VI::ILbr::KickOffDiscoveryOffset], 8);
}
// ILBR, branch 10
{
const auto* ilbr = &amswd.getILBr()[ilbrOffset(10)];
BOOST_CHECK_EQUAL(ilbr[VI::ILbr::OutletSegment], 23);
BOOST_CHECK_EQUAL(ilbr[VI::ILbr::NumBranchSegments], 1);
BOOST_CHECK_EQUAL(ilbr[VI::ILbr::FirstSegment], 29);
BOOST_CHECK_EQUAL(ilbr[VI::ILbr::LastSegment], 29);
BOOST_CHECK_EQUAL(ilbr[VI::ILbr::KickOffDiscoveryOffset], 9);
}
}
// The segments must appear in the following depth first search toe-to-heel
// order in ISEG[0]. We furthermore, go along kick-off branches before
// searching the main branch. Note that this order is *different* from
// ILBS/ILBR.
//
// 27 -- Branch ( 8)
// 28 -- Branch ( 9)
// 29 -- Branch (10)
// 23, 22, 21 -- Branch ( 5)
// 25 -- Branch ( 7)
// 24 -- Branch ( 6)
// 18, 17, 16, 15, 14, 13 -- Branch ( 3)
// 26, 20, 19 -- Branch ( 4)
// 12, 11, 10, 9, 8, 7 -- Branch ( 2)
// 6, 5, 4, 3, 2, 1 -- Branch ( 1)
//
BOOST_AUTO_TEST_CASE(Multilateral_ICD_Valve_ISEG_0)
{
const auto cse = SimulationCase { multilaterals_with_icd_valve() };
const auto& es = cse.es;
const auto& grid = cse.grid;
const auto& sched = cse.sched;
// Report Step 1: 2023-11-23 --> 2023-12-23
const auto rptStep = std::size_t {1};
const auto secs_elapsed = 30 * 86'400.0;
const auto ih = Opm::RestartIO::Helpers::
createInteHead(es, grid, sched, secs_elapsed,
rptStep, rptStep + 1, rptStep);
// No dynamic data needed for this test. We're checking the structure
// only.
const auto smry = Opm::SummaryState { Opm::TimeService::now() };
const auto xw = Opm::data::Wells {};
auto amswd = Opm::RestartIO::Helpers::AggregateMSWData {ih};
amswd.captureDeclaredMSWData(sched, rptStep, es.getUnits(),
ih, grid, smry, xw);
auto isegOffset = [&ih](const int ix)
{
return ih[VI::intehead::NISEGZ] * ix;
};
const auto expect = std::vector {
27, // Branch ( 8)
28, // Branch ( 9)
29, // Branch (10)
23, 22, 21, // Branch ( 5)
25, // Branch ( 7)
24, // Branch ( 6)
18, 17, 16, 15, 14, 13, // Branch ( 3)
26, 20, 19, // Branch ( 4)
12, 11, 10, 9, 8, 7, // Branch ( 2)
6, 5, 4, 3, 2, 1, // Branch ( 1)
};
const auto& iseg = amswd.getISeg();
for (auto i = 0*expect.size(); i < expect.size(); ++i) {
BOOST_CHECK_MESSAGE(iseg[isegOffset(i)] == expect[i],
"ISEG[0](" << i << ") == "
<< iseg[isegOffset(i)]
<< " differs from expected value "
<< expect[i]);
}
}
BOOST_AUTO_TEST_CASE(MSW_AICD)
{
const auto simCase = SimulationCase {first_sim("TEST_AGGREGATE_MSW.DATA")};