OilfieldToolsNet/opm-common
4
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Merge pull request #1559 from joakim-hove/rst-msw

Browse Source 
Initialize MSW wells from restart file
This commit is contained in:
Joakim Hove 2020-03-28 09:47:46 +01:00 committed by GitHub
commit 0025144079
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
14 changed files with 104 additions and 69 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

8
ExtraTests.cmake
View File

@ -66,10 +66,10 @@ opm_add_test( rst_spe1
LIBRARIES ${TEST_LIBS}
TEST_ARGS tests/SPE1CASE2.DATA tests/SPE1CASE2_RESTART.DATA )
# opm_add_test( rst_msw
# SOURCES tests/rst_test.cpp
# LIBRARIES ${TEST_LIBS}
# TEST_ARGS tests/MSW.DATA tests/MSW_RESTART.DATA )
opm_add_test( rst_msw
SOURCES tests/rst_test.cpp
LIBRARIES ${TEST_LIBS}
TEST_ARGS tests/MSW.DATA tests/MSW_RESTART.DATA )
# opm-tests dependent tests
if(HAVE_OPM_TESTS)

4
opm/io/eclipse/rst/connection.hpp
View File

@ -31,10 +31,10 @@ namespace RestartIO {
class Header;
struct RstConnection {
RstConnection(const ::Opm::UnitSystem& unit_system, const int* icon, const float* scon, const double *xcon);
RstConnection(const ::Opm::UnitSystem& unit_system, std::size_t rst_index, int nsconz, const int* icon, const float* scon, const double *xcon);
static double inverse_peaceman(double cf, double kh, double rw, double skin);
int insert_index;
std::size_t rst_index;
std::array<int,3> ijk;
Connection::State state;
int drain_sat_table;

2
opm/io/eclipse/rst/segment.hpp
View File

@ -30,7 +30,7 @@ class UnitSystem;
namespace RestartIO {
struct RstSegment {
RstSegment(const ::Opm::UnitSystem& unit_system, const int* iseg, const double * rseg);
RstSegment(const ::Opm::UnitSystem& unit_system, int segment_number, const int* iseg, const double * rseg);
int segment;
int outlet_segment;

4
opm/parser/eclipse/EclipseState/Schedule/Well/Connection.hpp
View File

@ -98,7 +98,7 @@ namespace RestartIO {
const double segDistEnd,
const bool defaultSatTabId);
Connection(const RestartIO::RstConnection& rst_connection, std::size_t insert_index, const EclipseGrid& grid, const FieldPropsManager& fp);
Connection(const RestartIO::RstConnection& rst_connection, const EclipseGrid& grid, const FieldPropsManager& fp);
static Connection serializeObject();
@ -125,6 +125,8 @@ namespace RestartIO {
void setState(State state);
void setComplnum(int compnum);
void scaleWellPi(double wellPi);
void updateSegmentRST(int segment_number_arg,
double center_depth_arg);
void updateSegment(int segment_number_arg,
double center_depth_arg,
std::size_t compseg_insert_index,

15
src/opm/io/eclipse/rst/connection.cpp
View File

@ -75,8 +75,8 @@ double RstConnection::inverse_peaceman(double cf, double kh, double rw, double s
using M = ::Opm::UnitSystem::measure;
RstConnection::RstConnection(const ::Opm::UnitSystem& unit_system, const int* icon, const float* scon, const double* xcon) :
insert_index( icon[VI::IConn::SeqIndex] - 1),
RstConnection::RstConnection(const ::Opm::UnitSystem& unit_system, std::size_t rst_index_, int nsconz, const int* icon, const float* scon, const double* xcon) :
rst_index( rst_index_),
ijk( {icon[VI::IConn::CellI] - 1, icon[VI::IConn::CellJ] - 1, icon[VI::IConn::CellK] - 1}),
state( from_int<Connection::State>(icon[VI::IConn::ConnStat])),
drain_sat_table( icon[VI::IConn::Drainage]),
@ -84,7 +84,7 @@ RstConnection::RstConnection(const ::Opm::UnitSystem& unit_system, const int* ic
completion( icon[VI::IConn::ComplNum]),
dir( from_int<Connection::Direction>(icon[VI::IConn::ConnDir])),
segment( icon[VI::IConn::Segment]),
cf_kind( from_float(scon[VI::SConn::CFInDeck])),
cf_kind( from_float(1.0)),
skin_factor( scon[VI::SConn::SkinFactor]),
cf( unit_system.to_si(M::transmissibility, scon[VI::SConn::ConnTrans])),
depth( unit_system.to_si(M::length, scon[VI::SConn::Depth])),
@ -103,12 +103,11 @@ RstConnection::RstConnection(const ::Opm::UnitSystem& unit_system, const int* ic
file. If the r0 value is given explicitly in the deck it is possible
to give a value which is not consistent with the Peaceman formula -
that value will be lost when loading back from a restart file.
insert_index: The insert index property is used internally to keep track
of the order the connections have been specified in the deck. It seems
that in some cases (MSW ?) eclipse only outputs 0 here.
*/
{ }
{
if (nsconz > VI::SConn::CFInDeck)
this->cf_kind = from_float(scon[VI::SConn::CFInDeck]);
}
}
}

4
src/opm/io/eclipse/rst/segment.cpp
View File

@ -50,8 +50,8 @@ double area_to_si(const UnitSystem& unit_system, double raw_value) {
}
RstSegment::RstSegment(const ::Opm::UnitSystem& unit_system, const int * iseg, const double * rseg) :
segment( iseg[VI::ISeg::SegNo]),
RstSegment::RstSegment(const ::Opm::UnitSystem& unit_system, int segment_number, const int * iseg, const double * rseg) :
segment( segment_number),
outlet_segment( iseg[VI::ISeg::OutSeg]),
branch( iseg[VI::ISeg::BranchNo]),
segment_type( from_ecl<Segment::SegmentType>(iseg[VI::ISeg::SegmentType])),

9
src/opm/io/eclipse/rst/well.cpp
View File

@ -109,7 +109,7 @@ RstWell::RstWell(const ::Opm::UnitSystem& unit_system,
std::size_t icon_offset = ic * header.niconz;
std::size_t scon_offset = ic * header.nsconz;
std::size_t xcon_offset = ic * header.nxconz;
this->connections.emplace_back( unit_system, icon + icon_offset, scon + scon_offset, xcon + xcon_offset);
this->connections.emplace_back( unit_system, ic, header.nsconz, icon + icon_offset, scon + scon_offset, xcon + xcon_offset);
}
}
@ -135,10 +135,11 @@ RstWell::RstWell(const ::Opm::UnitSystem& unit_system,
for (int is=0; is < header.nsegmx; is++) {
std::size_t iseg_offset = header.nisegz * (is + (this->msw_index - 1) * header.nsegmx);
std::size_t rseg_offset = header.nrsegz * (is + (this->msw_index - 1) * header.nsegmx);
auto segment_number = iseg[iseg_offset + VI::ISeg::SegNo];
if (segment_number != 0) {
auto other_segment_number = iseg[iseg_offset + VI::ISeg::SegNo];
if (other_segment_number != 0) {
auto segment_number = is + 1;
segment_map.insert({segment_number, this->segments.size()});
this->segments.emplace_back( unit_system, iseg.data() + iseg_offset, rseg.data() + rseg_offset);
this->segments.emplace_back( unit_system, segment_number, iseg.data() + iseg_offset, rseg.data() + rseg_offset);
}
}

17
src/opm/parser/eclipse/EclipseState/Schedule/MSW/Segment.cpp
View File

@ -47,6 +47,17 @@ static constexpr double invalid_value = -1.e100;
{
}
namespace {
double if_invalid_value(double rst_value) {
if (rst_value == 0)
return invalid_value;
return rst_value;
}
}
Segment::Segment(const RestartIO::RstSegment& rst_segment):
m_segment_number(rst_segment.segment),
@ -54,9 +65,9 @@ static constexpr double invalid_value = -1.e100;
m_outlet_segment(rst_segment.outlet_segment),
m_total_length( rst_segment.dist_bhp_ref ),
m_depth(rst_segment.bhp_ref_dz),
m_internal_diameter(rst_segment.diameter),
m_roughness(rst_segment.roughness),
m_cross_area(rst_segment.area),
m_internal_diameter(if_invalid_value(rst_segment.diameter)),
m_roughness(if_invalid_value(rst_segment.roughness)),
m_cross_area(if_invalid_value(rst_segment.area)),
m_volume(rst_segment.volume),
m_data_ready(true),
m_segment_type(rst_segment.segment_type)

93
src/opm/parser/eclipse/EclipseState/Schedule/Schedule.cpp
View File

@ -2958,7 +2958,7 @@ void Schedule::load_rst(const RestartIO::RstState& rst_state, const EclipseGrid&
std::unordered_map<int, Opm::Segment> segments;
for (const auto& rst_conn : rst_well.connections)
connections.emplace_back(rst_conn, connections.size(), grid, fp);
connections.emplace_back(rst_conn, grid, fp);
for (const auto& rst_segment : rst_well.segments) {
Opm::Segment segment(rst_segment);
@ -2968,15 +2968,9 @@ void Schedule::load_rst(const RestartIO::RstState& rst_state, const EclipseGrid&
for (auto& connection : connections) {
int segment_id = connection.segment();
if (segment_id > 0) {
std::size_t compsegs_insert_index = 0;
double segment_start = 0;
double segment_end = 0;
const auto& segment = segments.at(segment_id);
connection.updateSegment(segment.segmentNumber(),
segment.depth(),
compsegs_insert_index,
segment_start,
segment_end);
connection.updateSegmentRST(segment.segmentNumber(),
segment.depth());
}
}
@ -2987,9 +2981,19 @@ void Schedule::load_rst(const RestartIO::RstState& rst_state, const EclipseGrid&
if (!segments.empty()) {
std::vector<Segment> segments_list;
/*
The ordering of the segments in the WellSegments structure seems a
bit random; in some parts of the code the segment_number seems to
be treated like a random integer ID, whereas in other parts it
seems to be treated like a running index. Here the segments in
WellSegments are sorted according to the segment number - observe
that this is somewhat important because the first top segment is
treated differently from the other segment.
*/
for (const auto& segment_pair : segments)
segments_list.push_back( std::move(segment_pair.second) );
std::sort( segments_list.begin(), segments_list.end(),[](const Segment& seg1, const Segment& seg2) { return seg1.segmentNumber() < seg2.segmentNumber(); } );
auto comp_pressure_drop = WellSegments::CompPressureDrop::HFA;
std::shared_ptr<Opm::WellSegments> well_segments = std::make_shared<Opm::WellSegments>(comp_pressure_drop, segments_list);
well.updateSegments( std::move(well_segments) );
@ -3069,6 +3073,13 @@ std::string well_msg(const std::string& well, const std::string& msg) {
return "Well: " + well + " " + msg;
}
std::string well_segment_msg(const std::string& well, int segment_number, const std::string& msg) {
return "Well: " + well + " Segment: " + std::to_string(segment_number) + " " + msg;
}
std::string well_connection_msg(const std::string& well, const Connection& conn, const std::string& msg) {
return "Well: " + well + " Connection: " + std::to_string(conn.getI()) + ", " + std::to_string(conn.getJ()) + ", " + std::to_string(conn.getK()) + " " + msg;
}
}
@ -3102,26 +3113,25 @@ bool Schedule::cmp(const Schedule& sched1, const Schedule& sched2, std::size_t r
for (std::size_t icon = 0; icon < connections1.size(); icon++) {
const auto& conn1 = connections1[icon];
const auto& conn2 = connections2[icon];
well_count += not_equal( conn1.getI(), conn2.getI(), well_msg(well1.name(), "Connection: I"));
well_count += not_equal( conn1.getI() , conn2.getI() , well_msg(well1.name(), "Connection: I"));
well_count += not_equal( conn1.getJ() , conn2.getJ() , well_msg(well1.name(), "Connection: J"));
well_count += not_equal( conn1.getK() , conn2.getK() , well_msg(well1.name(), "Connection: K"));
well_count += not_equal( conn1.state() , conn2.state(), well_msg(well1.name(), "Connection: State"));
well_count += not_equal( conn1.dir() , conn2.dir(), well_msg(well1.name(), "Connection: dir"));
well_count += not_equal( conn1.complnum() , conn2.complnum(), well_msg(well1.name(), "connection: complnum"));
well_count += not_equal( conn1.segment() , conn2.segment(), well_msg(well1.name(), "Connection: segment"));
well_count += not_equal( conn1.kind() , conn2.kind(), well_msg(well1.name(), "Connection: CFKind"));
well_count += not_equal( conn1.sort_value(), conn2.sort_value(), well_msg(well1.name(), "Connection: sort_value"));
well_count += not_equal( conn1.getI(), conn2.getI(), well_connection_msg(well1.name(), conn1, "I"));
well_count += not_equal( conn1.getJ() , conn2.getJ() , well_connection_msg(well1.name(), conn1, "J"));
well_count += not_equal( conn1.getK() , conn2.getK() , well_connection_msg(well1.name(), conn1, "K"));
well_count += not_equal( conn1.state() , conn2.state(), well_connection_msg(well1.name(), conn1, "State"));
well_count += not_equal( conn1.dir() , conn2.dir(), well_connection_msg(well1.name(), conn1, "dir"));
well_count += not_equal( conn1.complnum() , conn2.complnum(), well_connection_msg(well1.name(), conn1, "complnum"));
well_count += not_equal( conn1.segment() , conn2.segment(), well_connection_msg(well1.name(), conn1, "segment"));
well_count += not_equal( conn1.kind() , conn2.kind(), well_connection_msg(well1.name(), conn1, "CFKind"));
well_count += not_equal( conn1.sort_value(), conn2.sort_value(), well_connection_msg(well1.name(), conn1, "sort_value"));
well_count += not_equal( conn1.CF(), conn2.CF(), well_msg(well1.name(), "Connection: CF"));
well_count += not_equal( conn1.Kh(), conn2.Kh(), well_msg(well1.name(), "Connection: Kh"));
well_count += not_equal( conn1.rw(), conn2.rw(), well_msg(well1.name(), "Connection: rw"));
well_count += not_equal( conn1.depth(), conn2.depth(), well_msg(well1.name(), "Connection: depth"));
well_count += not_equal( conn1.CF(), conn2.CF(), well_connection_msg(well1.name(), conn1, "CF"));
well_count += not_equal( conn1.Kh(), conn2.Kh(), well_connection_msg(well1.name(), conn1, "Kh"));
well_count += not_equal( conn1.rw(), conn2.rw(), well_connection_msg(well1.name(), conn1, "rw"));
well_count += not_equal( conn1.depth(), conn2.depth(), well_connection_msg(well1.name(), conn1, "depth"));
well_count += not_equal( conn1.r0(), conn2.r0(), well_msg(well1.name(), "Connection: r0"));
well_count += not_equal( conn1.skinFactor(), conn2.skinFactor(), well_msg(well1.name(), "Connection: skinFactor"));
well_count += not_equal( conn1.wellPi(), conn2.wellPi(), well_msg(well1.name(), "Connection: PI"));
//well_count += not_equal( conn1.r0(), conn2.r0(), well_connection_msg(well1.name(), conn1, "r0"));
well_count += not_equal( conn1.skinFactor(), conn2.skinFactor(), well_connection_msg(well1.name(), conn1, "skinFactor"));
well_count += not_equal( conn1.wellPi(), conn2.wellPi(), well_connection_msg(well1.name(), conn1, "PI"));
}
}
@ -3138,15 +3148,16 @@ bool Schedule::cmp(const Schedule& sched1, const Schedule& sched2, std::size_t r
for (std::size_t iseg=0; iseg < segments1.size(); iseg++) {
const auto& segment1 = segments1[iseg];
const auto& segment2 = segments2[iseg];
well_count += not_equal(segment1.segmentNumber(), segment2.segmentNumber(), well_msg(well1.name(), "Segment: segmentNumber"));
well_count += not_equal(segment1.branchNumber(), segment2.segmentNumber(), well_msg(well1.name(), "Segment: segmentNumber"));
well_count += not_equal(segment1.outletSegment(), segment2.outletSegment(), well_msg(well1.name(), "Segments: outletSegment"));
well_count += not_equal(segment1.totalLength(), segment2.totalLength(), well_msg(well1.name(), "Segments: totalLength"));
well_count += not_equal(segment1.depth(), segment2.depth(), well_msg(well1.name(), "Segments: depth"));
well_count += not_equal(segment1.internalDiameter(), segment2.internalDiameter(), well_msg(well1.name(), "Segments: internalDiameter"));
well_count += not_equal(segment1.roughness(), segment2.roughness(), well_msg(well1.name(), "Segments: roughness"));
well_count += not_equal(segment1.crossArea(), segment2.crossArea(), well_msg(well1.name(), "Segments: crossArea"));
well_count += not_equal(segment1.volume(), segment2.volume(), well_msg(well1.name(), "Segments: volume"));
//const auto& segment2 = segments2.getFromSegmentNumber(segment1.segmentNumber());
well_count += not_equal(segment1.segmentNumber(), segment2.segmentNumber(), well_segment_msg(well1.name(), segment1.segmentNumber(), "segmentNumber"));
well_count += not_equal(segment1.branchNumber(), segment2.branchNumber(), well_segment_msg(well1.name(), segment1.segmentNumber(), "branchNumber"));
well_count += not_equal(segment1.outletSegment(), segment2.outletSegment(), well_segment_msg(well1.name(), segment1.segmentNumber(), "outletSegment"));
well_count += not_equal(segment1.totalLength(), segment2.totalLength(), well_segment_msg(well1.name(), segment1.segmentNumber(), "totalLength"));
well_count += not_equal(segment1.depth(), segment2.depth(), well_segment_msg(well1.name(), segment1.segmentNumber(), "depth"));
well_count += not_equal(segment1.internalDiameter(), segment2.internalDiameter(), well_segment_msg(well1.name(), segment1.segmentNumber(), "internalDiameter"));
well_count += not_equal(segment1.roughness(), segment2.roughness(), well_segment_msg(well1.name(), segment1.segmentNumber(), "roughness"));
well_count += not_equal(segment1.crossArea(), segment2.crossArea(), well_segment_msg(well1.name(), segment1.segmentNumber(), "crossArea"));
well_count += not_equal(segment1.volume(), segment2.volume(), well_segment_msg(well1.name(), segment1.segmentNumber(), "volume"));
}
}
@ -3162,14 +3173,16 @@ bool Schedule::cmp(const Schedule& sched1, const Schedule& sched2, std::size_t r
well_count += not_equal(prod1.ResVRate, prod2.ResVRate, well_msg(well1.name(), "Prod: ResVRate"));
well_count += not_equal(prod1.BHPTarget, prod2.BHPTarget, well_msg(well1.name(), "Prod: BHPTarget"));
well_count += not_equal(prod1.THPTarget, prod2.THPTarget, well_msg(well1.name(), "Prod: THPTarget"));
well_count += not_equal(prod1.bhp_hist_limit, prod2.bhp_hist_limit, well_msg(well1.name(), "Prod: bhp_hist_limit"));
well_count += not_equal(prod1.thp_hist_limit, prod2.thp_hist_limit, well_msg(well1.name(), "Prod: thp_hist_limit"));
well_count += not_equal(prod1.BHPH, prod2.BHPH, well_msg(well1.name(), "Prod: BHPH"));
well_count += not_equal(prod1.THPH, prod2.THPH, well_msg(well1.name(), "Prod: THPH"));
well_count += not_equal(prod1.predictionMode, prod2.predictionMode, well_msg(well1.name(), "Prod: predictionMode"));
if (!prod1.predictionMode) {
well_count += not_equal(prod1.bhp_hist_limit, prod2.bhp_hist_limit, well_msg(well1.name(), "Prod: bhp_hist_limit"));
well_count += not_equal(prod1.thp_hist_limit, prod2.thp_hist_limit, well_msg(well1.name(), "Prod: thp_hist_limit"));
well_count += not_equal(prod1.BHPH, prod2.BHPH, well_msg(well1.name(), "Prod: BHPH"));
well_count += not_equal(prod1.THPH, prod2.THPH, well_msg(well1.name(), "Prod: THPH"));
}
well_count += not_equal(prod1.VFPTableNumber, prod2.VFPTableNumber, well_msg(well1.name(), "Prod: VFPTableNumber"));
well_count += not_equal(prod1.ALQValue, prod2.ALQValue, well_msg(well1.name(), "Prod: ALQValue"));
well_count += not_equal(prod1.productionControls(), prod2.productionControls(), well_msg(well1.name(), "Prod: productionControls"));
well_count += not_equal(prod1.predictionMode, prod2.predictionMode, well_msg(well1.name(), "Prod: predictionMode"));
well_count += not_equal(prod1.controlMode, prod2.controlMode, well_msg(well1.name(), "Prod: controlMode"));
well_count += not_equal(prod1.whistctl_cmode, prod2.whistctl_cmode, well_msg(well1.name(), "Prod: whistctl_cmode"));
}

10
src/opm/parser/eclipse/EclipseState/Schedule/Well/Connection.cpp
View File

@ -82,7 +82,7 @@ constexpr bool defaultSatTabId = true;
constexpr double def_wellPi = 1.0;
}
Connection::Connection(const RestartIO::RstConnection& rst_connection, std::size_t insert_index, const EclipseGrid& grid, const FieldPropsManager& fp) :
Connection::Connection(const RestartIO::RstConnection& rst_connection, const EclipseGrid& grid, const FieldPropsManager& fp) :
direction(rst_connection.dir),
center_depth(rst_connection.depth),
open_state(rst_connection.state),
@ -96,7 +96,7 @@ Connection::Connection(const RestartIO::RstConnection& rst_connection, std::size
ijk(rst_connection.ijk),
m_ctfkind(rst_connection.cf_kind),
m_global_index(grid.getGlobalIndex(this->ijk[0], this->ijk[1], this->ijk[2])),
m_sort_value(insert_index),
m_sort_value(rst_connection.rst_index),
m_segDistStart(rst_connection.segdist_start),
m_segDistEnd(rst_connection.segdist_end),
m_defaultSatTabId(defaultSatTabId),
@ -250,6 +250,12 @@ Connection::Connection(const RestartIO::RstConnection& rst_connection, std::size
this->m_segDistEnd = end;
}
void Connection::updateSegmentRST(int segment_number_arg,
double center_depth_arg) {
this->segment_number = segment_number_arg;
this->center_depth = center_depth_arg;
}
int Connection::segment() const {
return this->segment_number;
}

2
tests/parser/ScheduleRestartTests.cpp
View File

@ -48,7 +48,7 @@ void compare_connections(const RestartIO::RstConnection& rst_conn, const Connect
BOOST_CHECK_EQUAL(rst_conn.ijk[2], sched_conn.getK());
BOOST_CHECK_EQUAL(rst_conn.segment, sched_conn.segment());
BOOST_CHECK_EQUAL(rst_conn.insert_index, static_cast<int>(sched_conn.sort_value()));
BOOST_CHECK_EQUAL(rst_conn.rst_index, static_cast<int>(sched_conn.sort_value()));
BOOST_CHECK(rst_conn.state == sched_conn.state());
BOOST_CHECK(rst_conn.dir == sched_conn.dir());
BOOST_CHECK_CLOSE( rst_conn.cf, sched_conn.CF() , 1e-6);

BIN
tests/restart/MSW.UNRST
View File

Binary file not shown.

3
tests/restart/README Normal file
View File

@ -0,0 +1,3 @@
The restart file 'MSW.UNRST' contains report step from the reference solution in
opm-tests/msw_3d_hfa. If the reference solution in opm-tests is updated the
restart file 'MSW.UNRST' must be updated manually :-(

2
tests/test_AggregateMSWData.cpp
View File

@ -845,7 +845,7 @@ BOOST_AUTO_TEST_CASE(MSW_RST) {
);
const auto& iseg = amswd.getISeg();
const auto& rseg = amswd.getRSeg();
auto segment = Opm::RestartIO::RstSegment(simCase.es.getUnits(), iseg.data(), rseg.data());
auto segment = Opm::RestartIO::RstSegment(simCase.es.getUnits(), 1, iseg.data(), rseg.data());
}