diff --git a/msim/include/opm/msim/msim.hpp b/msim/include/opm/msim/msim.hpp index e22783730..27e97a500 100644 --- a/msim/include/opm/msim/msim.hpp +++ b/msim/include/opm/msim/msim.hpp @@ -12,6 +12,7 @@ #include #include #include +#include #include #include @@ -29,10 +30,6 @@ class SummaryState; class UDQState; class WellTestState; -namespace Action { -class State; -} - class msim { public: @@ -46,13 +43,16 @@ public: void well_rate(const std::string& well, data::Rates::opt rate, std::function func); void solution(const std::string& field, std::function func); void run(Schedule& schedule, EclipseIO& io, bool report_only); - void post_step(Schedule& schedule, Action::State& action_state, SummaryState& st, data::Solution& sol, data::Wells& well_data, data::GroupAndNetworkValues& group_nwrk_data, size_t report_step, const time_point& sim_time); + void post_step(Schedule& schedule, SummaryState& st, data::Solution& sol, data::Wells& well_data, data::GroupAndNetworkValues& group_nwrk_data, size_t report_step, const time_point& sim_time); + + Action::State action_state; + private: - void run_step(const Schedule& schedule, Action::State& action_state, WellTestState& wtest_state, SummaryState& st, UDQState& udq_state, data::Solution& sol, data::Wells& well_data, data::GroupAndNetworkValues& group_nwrk_data, size_t report_step, EclipseIO& io) const; - void run_step(const Schedule& schedule, Action::State& action_state, WellTestState& wtest_state, SummaryState& st, UDQState& udq_state, data::Solution& sol, data::Wells& well_data, data::GroupAndNetworkValues& group_nwrk_data, size_t report_step, double dt, EclipseIO& io) const; - void output(Action::State& action_state, WellTestState& wtest_state, SummaryState& st, const UDQState& udq_state, size_t report_step, bool substep, double seconds_elapsed, const data::Solution& sol, const data::Wells& well_data, const data::GroupAndNetworkValues& group_data, EclipseIO& io) const; - void simulate(const Schedule& schedule, WellTestState& wtest_state, const SummaryState& st, data::Solution& sol, data::Wells& well_data, data::GroupAndNetworkValues& group_nwrk_data, size_t report_step, double seconds_elapsed, double time_step) const; + void run_step(const Schedule& schedule, WellTestState& wtest_state, SummaryState& st, UDQState& udq_state, data::Solution& sol, data::Wells& well_data, data::GroupAndNetworkValues& group_nwrk_data, size_t report_step, EclipseIO& io) const; + void run_step(const Schedule& schedule, WellTestState& wtest_state, SummaryState& st, UDQState& udq_state, data::Solution& sol, data::Wells& well_data, data::GroupAndNetworkValues& group_nwrk_data, size_t report_step, double dt, EclipseIO& io) const; + void output(WellTestState& wtest_state, SummaryState& st, const UDQState& udq_state, size_t report_step, bool substep, double seconds_elapsed, const data::Solution& sol, const data::Wells& well_data, const data::GroupAndNetworkValues& group_data, EclipseIO& io) const; + void simulate(const Schedule& schedule, const SummaryState& st, data::Solution& sol, data::Wells& well_data, data::GroupAndNetworkValues& group_nwrk_data, size_t report_step, double seconds_elapsed, double time_step) const; EclipseState state; std::map>> well_rates; diff --git a/msim/src/msim.cpp b/msim/src/msim.cpp index 681c456f5..769808932 100644 --- a/msim/src/msim.cpp +++ b/msim/src/msim.cpp @@ -31,7 +31,6 @@ #include #include -#include #include #include #include @@ -52,7 +51,6 @@ void msim::run(Schedule& schedule, EclipseIO& io, bool report_only) { SummaryState st(TimeService::from_time_t(schedule.getStartTime())); UDQState udq_state(schedule.getUDQConfig(0).params().undefinedValue()); WellTestState wtest_state; - Action::State action_state; Python python; io.writeInitial(); @@ -60,13 +58,13 @@ void msim::run(Schedule& schedule, EclipseIO& io, bool report_only) { data::Wells well_data; data::GroupAndNetworkValues group_nwrk_data; if (report_only) - run_step(schedule, action_state, wtest_state, st, udq_state, sol, well_data, group_nwrk_data, report_step, io); + run_step(schedule, wtest_state, st, udq_state, sol, well_data, group_nwrk_data, report_step, io); else { double time_step = std::min(week, 0.5*schedule.stepLength(report_step - 1)); - run_step(schedule, action_state, wtest_state, st, udq_state, sol, well_data, group_nwrk_data, report_step, time_step, io); + run_step(schedule, wtest_state, st, udq_state, sol, well_data, group_nwrk_data, report_step, time_step, io); } auto sim_time = TimeService::from_time_t( schedule.simTime(report_step) ); - post_step(schedule, action_state, st, sol, well_data, group_nwrk_data, report_step, sim_time); + post_step(schedule, st, sol, well_data, group_nwrk_data, report_step, sim_time); const auto& exit_status = schedule.exitStatus(); if (exit_status.has_value()) return; @@ -78,14 +76,14 @@ UDAValue msim::uda_val() { } -void msim::post_step(Schedule& schedule, Action::State& action_state, SummaryState& st, data::Solution& /* sol */, data::Wells& /* well_data */, data::GroupAndNetworkValues& /* grp_nwrk_data */, size_t report_step, const time_point& sim_time) { +void msim::post_step(Schedule& schedule, SummaryState& st, data::Solution& /* sol */, data::Wells& /* well_data */, data::GroupAndNetworkValues& /* grp_nwrk_data */, size_t report_step, const time_point& sim_time) { const auto& actions = schedule[report_step].actions.get(); if (actions.empty()) return; Action::Context context( st , schedule[report_step].wlist_manager.get()); - for (const auto& action : actions.pending(action_state, std::chrono::system_clock::to_time_t(sim_time))) { + for (const auto& action : actions.pending(this->action_state, std::chrono::system_clock::to_time_t(sim_time))) { auto result = action->eval(context); if (result) schedule.applyAction(report_step, sim_time, *action, result, {}); @@ -97,12 +95,12 @@ void msim::post_step(Schedule& schedule, Action::State& action_state, SummarySta -void msim::run_step(const Schedule& schedule, Action::State& action_state, WellTestState& wtest_state, SummaryState& st, UDQState& udq_state, data::Solution& sol, data::Wells& well_data, data::GroupAndNetworkValues& grp_nwrk_data, size_t report_step, EclipseIO& io) const { - this->run_step(schedule, action_state, wtest_state, st, udq_state, sol, well_data, grp_nwrk_data, report_step, schedule.stepLength(report_step - 1), io); +void msim::run_step(const Schedule& schedule, WellTestState& wtest_state, SummaryState& st, UDQState& udq_state, data::Solution& sol, data::Wells& well_data, data::GroupAndNetworkValues& grp_nwrk_data, size_t report_step, EclipseIO& io) const { + this->run_step(schedule, wtest_state, st, udq_state, sol, well_data, grp_nwrk_data, report_step, schedule.stepLength(report_step - 1), io); } -void msim::run_step(const Schedule& schedule, Action::State& action_state, WellTestState& wtest_state, SummaryState& st, UDQState& udq_state, data::Solution& sol, data::Wells& well_data, data::GroupAndNetworkValues& group_nwrk_data, size_t report_step, double dt, EclipseIO& io) const { +void msim::run_step(const Schedule& schedule, WellTestState& wtest_state, SummaryState& st, UDQState& udq_state, data::Solution& sol, data::Wells& well_data, data::GroupAndNetworkValues& group_nwrk_data, size_t report_step, double dt, EclipseIO& io) const { double start_time = schedule.seconds(report_step - 1); double end_time = schedule.seconds(report_step); double seconds_elapsed = start_time; @@ -112,7 +110,7 @@ void msim::run_step(const Schedule& schedule, Action::State& action_state, WellT if ((seconds_elapsed + time_step) > end_time) time_step = end_time - seconds_elapsed; - this->simulate(schedule, wtest_state, st, sol, well_data, group_nwrk_data, report_step, seconds_elapsed, time_step); + this->simulate(schedule, st, sol, well_data, group_nwrk_data, report_step, seconds_elapsed, time_step); seconds_elapsed += time_step; @@ -128,8 +126,7 @@ void msim::run_step(const Schedule& schedule, Action::State& action_state, WellT schedule.getUDQConfig( report_step ).eval(report_step, schedule.wellMatcher(report_step), st, udq_state); - this->output(action_state, - wtest_state, + this->output(wtest_state, st, udq_state, report_step, @@ -144,9 +141,9 @@ void msim::run_step(const Schedule& schedule, Action::State& action_state, WellT -void msim::output(Action::State& action_state, WellTestState& wtest_state, SummaryState& st, const UDQState& udq_state, size_t report_step, bool substep, double seconds_elapsed, const data::Solution& sol, const data::Wells& well_data, const data::GroupAndNetworkValues& group_nwrk_data, EclipseIO& io) const { +void msim::output(WellTestState& wtest_state, SummaryState& st, const UDQState& udq_state, size_t report_step, bool substep, double seconds_elapsed, const data::Solution& sol, const data::Wells& well_data, const data::GroupAndNetworkValues& group_nwrk_data, EclipseIO& io) const { RestartValue value(sol, well_data, group_nwrk_data, {}); - io.writeTimeStep(action_state, + io.writeTimeStep(this->action_state, wtest_state, st, udq_state, @@ -157,7 +154,7 @@ void msim::output(Action::State& action_state, WellTestState& wtest_state, Summa } -void msim::simulate(const Schedule& schedule, WellTestState&, const SummaryState& st, data::Solution& sol, data::Wells& well_data, data::GroupAndNetworkValues& /* group_nwrk_data */, size_t report_step, double seconds_elapsed, double time_step) const { +void msim::simulate(const Schedule& schedule, const SummaryState& st, data::Solution& sol, data::Wells& well_data, data::GroupAndNetworkValues& /* group_nwrk_data */, size_t report_step, double seconds_elapsed, double time_step) const { for (const auto& sol_pair : this->solutions) { auto func = sol_pair.second; func(this->state, schedule, sol, report_step, seconds_elapsed + time_step); diff --git a/opm/parser/eclipse/EclipseState/Schedule/CompletedCells.hpp b/opm/parser/eclipse/EclipseState/Schedule/CompletedCells.hpp index 977e1755d..d3e790fbb 100644 --- a/opm/parser/eclipse/EclipseState/Schedule/CompletedCells.hpp +++ b/opm/parser/eclipse/EclipseState/Schedule/CompletedCells.hpp @@ -37,15 +37,41 @@ public: double permx; double permy; double permz; - //int satnum; + int satnum; + int pvtnum; + double ntg; bool operator==(const Props& other) const{ return this->active_index == other.active_index && this->permx == other.permx && this->permy == other.permy && - this->permz == other.permz; + this->permz == other.permz && + this->satnum == other.satnum && + this->pvtnum == other.pvtnum && + this->ntg == other.ntg; } + template + void serializeOp(Serializer& serializer) + { + serializer(this->permx); + serializer(this->permy); + serializer(this->permz); + serializer(this->satnum); + serializer(this->pvtnum); + serializer(this->ntg); + } + + static Props serializeObject(){ + Props props; + props.permx = 10.0; + props.permy = 78.0; + props.permz = 45.4; + props.satnum = 3; + props.pvtnum = 5; + props.ntg = 45.1; + return props; + } }; std::optional props; diff --git a/opm/parser/eclipse/EclipseState/Schedule/Schedule.hpp b/opm/parser/eclipse/EclipseState/Schedule/Schedule.hpp index 6cabb3ebe..74244daf5 100644 --- a/opm/parser/eclipse/EclipseState/Schedule/Schedule.hpp +++ b/opm/parser/eclipse/EclipseState/Schedule/Schedule.hpp @@ -506,7 +506,6 @@ namespace Opm ErrorGuard& errors, const ScheduleGrid& grid, const std::unordered_map * target_wellpi, - const FieldPropsManager* fp, const std::string& prefix); void addACTIONX(const Action::ActionX& action); void addGroupToGroup( const std::string& parent_group, const std::string& child_group); @@ -521,7 +520,6 @@ namespace Opm const DeckKeyword& keyword, const ParseContext& parseContext, ErrorGuard& errors, const ScheduleGrid& grid, - const FieldPropsManager* fp, const std::vector& matching_wells, bool runtime, SimulatorUpdate * sim_update, @@ -544,7 +542,6 @@ namespace Opm SimulatorUpdate * sim_update; const std::unordered_map * target_wellpi; const ScheduleGrid& grid; - const FieldPropsManager* fp_ptr; HandlerContext(const ScheduleBlock& block_, const DeckKeyword& keyword_, @@ -562,7 +559,6 @@ namespace Opm , sim_update(sim_update_) , target_wellpi(target_wellpi_) , grid(grid_) - , fp_ptr(nullptr) {} void affected_well(const std::string& well_name) { diff --git a/opm/parser/eclipse/EclipseState/Schedule/Well/Well.hpp b/opm/parser/eclipse/EclipseState/Schedule/Well/Well.hpp index 416f52425..239a2d581 100644 --- a/opm/parser/eclipse/EclipseState/Schedule/Well/Well.hpp +++ b/opm/parser/eclipse/EclipseState/Schedule/Well/Well.hpp @@ -586,7 +586,7 @@ public: bool updateDrainageRadius(double drainage_radius); void updateSegments(std::shared_ptr segments_arg); bool updateConnections(std::shared_ptr connections, bool force); - bool updateConnections(std::shared_ptr connections, const ScheduleGrid& grid, const std::vector& pvtnum); + bool updateConnections(std::shared_ptr connections, const ScheduleGrid& grid); bool updateStatus(Status status); bool updateGroup(const std::string& group); bool updateWellGuideRate(bool available, double guide_rate, GuideRateTarget guide_phase, double scale_factor); diff --git a/opm/parser/eclipse/EclipseState/Schedule/Well/WellConnections.hpp b/opm/parser/eclipse/EclipseState/Schedule/Well/WellConnections.hpp index f007ad176..6ea10d0b5 100644 --- a/opm/parser/eclipse/EclipseState/Schedule/Well/WellConnections.hpp +++ b/opm/parser/eclipse/EclipseState/Schedule/Well/WellConnections.hpp @@ -80,7 +80,7 @@ namespace Opm { const Connection::CTFKind ctf_kind = Connection::CTFKind::DeckValue, const std::size_t seqIndex = 0, const bool defaultSatTabId = true); - void loadCOMPDAT(const DeckRecord& record, const ScheduleGrid& grid, const FieldPropsManager& field_properties, const std::string& wname, const KeywordLocation& location); + void loadCOMPDAT(const DeckRecord& record, const ScheduleGrid& grid, const std::string& wname, const KeywordLocation& location); using const_iterator = std::vector< Connection >::const_iterator; @@ -169,16 +169,6 @@ namespace Opm { const std::size_t seqIndex = 0, const bool defaultSatTabId = true); - void loadCOMPDAT(const DeckRecord& record, - const ScheduleGrid& grid, - const std::vector& satnum_data, - const std::vector* permx, - const std::vector* permy, - const std::vector* permz, - const std::vector& ntg, - const std::string& wname, - const KeywordLocation& location); - size_t findClosestConnection(int oi, int oj, double oz, size_t start_pos); void orderTRACK(); void orderMSW(); diff --git a/src/opm/parser/eclipse/EclipseState/Schedule/KeywordHandlers.cpp b/src/opm/parser/eclipse/EclipseState/Schedule/KeywordHandlers.cpp index c864a11ae..63b486d2d 100644 --- a/src/opm/parser/eclipse/EclipseState/Schedule/KeywordHandlers.cpp +++ b/src/opm/parser/eclipse/EclipseState/Schedule/KeywordHandlers.cpp @@ -155,8 +155,8 @@ namespace { for (const auto& name : wellnames) { auto well2 = this->snapshots.back().wells.get(name); auto connections = std::shared_ptr( new WellConnections( well2.getConnections())); - connections->loadCOMPDAT(record, handlerContext.grid, *handlerContext.fp_ptr, name, handlerContext.keyword.location()); - if (well2.updateConnections(connections, handlerContext.grid, handlerContext.fp_ptr->get_int("PVTNUM"))) { + connections->loadCOMPDAT(record, handlerContext.grid, name, handlerContext.keyword.location()); + if (well2.updateConnections(connections, handlerContext.grid)) { this->snapshots.back().wells.update( well2 ); wells.insert( name ); } diff --git a/src/opm/parser/eclipse/EclipseState/Schedule/Schedule.cpp b/src/opm/parser/eclipse/EclipseState/Schedule/Schedule.cpp index 7e5650d39..1412f7dc3 100644 --- a/src/opm/parser/eclipse/EclipseState/Schedule/Schedule.cpp +++ b/src/opm/parser/eclipse/EclipseState/Schedule/Schedule.cpp @@ -167,13 +167,13 @@ namespace Opm { if (rst) { auto restart_step = this->m_static.rst_info.report_step; - this->iterateScheduleSection( 0, restart_step, parseContext, errors, grid, nullptr, &fp, ""); + this->iterateScheduleSection( 0, restart_step, parseContext, errors, grid, nullptr, ""); this->load_rst(*rst, grid, fp); if (! this->restart_output.writeRestartFile(restart_step)) this->restart_output.addRestartOutput(restart_step); - this->iterateScheduleSection( restart_step, this->m_sched_deck.size(), parseContext, errors, grid, nullptr, &fp, ""); + this->iterateScheduleSection( restart_step, this->m_sched_deck.size(), parseContext, errors, grid, nullptr, ""); } else { - this->iterateScheduleSection( 0, this->m_sched_deck.size(), parseContext, errors, grid, nullptr, &fp, ""); + this->iterateScheduleSection( 0, this->m_sched_deck.size(), parseContext, errors, grid, nullptr, ""); } //m_grid = std::make_shared(grid, gridWrapper.getHitKeys()); @@ -297,7 +297,6 @@ Schedule::Schedule(const Deck& deck, const EclipseState& es, const std::optional const ParseContext& parseContext, ErrorGuard& errors, const ScheduleGrid& grid, - const FieldPropsManager* fp, const std::vector& matching_wells, bool actionx_mode, SimulatorUpdate * sim_update, @@ -314,9 +313,6 @@ Schedule::Schedule(const Deck& deck, const EclipseState& es, const std::optional The grid and fieldProps members create problems for reiterating the Schedule section. We therefor single them out very clearly here. */ - if (require_grid.count(keyword.name()) > 0) { - handlerContext.fp_ptr = fp; - } if (handleNormalKeyword(handlerContext, parseContext, errors)) return; @@ -392,7 +388,6 @@ void Schedule::iterateScheduleSection(std::size_t load_start, std::size_t load_e ErrorGuard& errors, const ScheduleGrid& grid, const std::unordered_map * target_wellpi, - const FieldPropsManager* fp, const std::string& prefix) { std::vector > rftProperties; @@ -507,7 +502,6 @@ void Schedule::iterateScheduleSection(std::size_t load_start, std::size_t load_e parseContext, errors, grid, - fp, {}, false, nullptr, @@ -1269,7 +1263,6 @@ void Schedule::iterateScheduleSection(std::size_t load_start, std::size_t load_e parseContext, errors, grid, - nullptr, result.wells(), true, &sim_update, @@ -1283,7 +1276,7 @@ void Schedule::iterateScheduleSection(std::size_t load_start, std::size_t load_e } if (reportStep < this->m_sched_deck.size() - 1) - iterateScheduleSection(reportStep + 1, this->m_sched_deck.size(), parseContext, errors, grid, &target_wellpi, nullptr, prefix); + iterateScheduleSection(reportStep + 1, this->m_sched_deck.size(), parseContext, errors, grid, &target_wellpi, prefix); OpmLog::info("\\----------------------------------------------------------------------"); return sim_update; @@ -1447,7 +1440,7 @@ namespace { rst_well.ij[0], rst_well.ij[1], rst_connections); - well.updateConnections( std::make_shared( std::move(connections) ), grid, fp.get_int("PVTNUM")); + well.updateConnections( std::make_shared( std::move(connections) ), grid); } else { std::unordered_map rst_segments; for (const auto& rst_segment : rst_well.segments) { @@ -1456,7 +1449,7 @@ namespace { } auto [connections, segments] = Compsegs::rstUpdate(rst_well, rst_connections, rst_segments); - well.updateConnections( std::make_shared(std::move(connections)), grid, fp.get_int("PVTNUM")); + well.updateConnections( std::make_shared(std::move(connections)), grid); well.updateSegments( std::make_shared(std::move(segments) )); } diff --git a/src/opm/parser/eclipse/EclipseState/Schedule/ScheduleGrid.cpp b/src/opm/parser/eclipse/EclipseState/Schedule/ScheduleGrid.cpp index f6158fe3d..eaed3d902 100644 --- a/src/opm/parser/eclipse/EclipseState/Schedule/ScheduleGrid.cpp +++ b/src/opm/parser/eclipse/EclipseState/Schedule/ScheduleGrid.cpp @@ -41,6 +41,13 @@ namespace { else throw std::logic_error(fmt::format("FieldPropsManager is missing keyword '{}'", kw)); } + + double try_get_ntg_value(const Opm::FieldPropsManager& fp, const std::string& kw, std::size_t active_index){ + if (fp.has_double(kw)) + return fp.try_get(kw)->at(active_index); + else + return 1.0; + } } const Opm::CompletedCells::Cell& Opm::ScheduleGrid::get_cell(std::size_t i, std::size_t j, std::size_t k) const { @@ -56,7 +63,9 @@ const Opm::CompletedCells::Cell& Opm::ScheduleGrid::get_cell(std::size_t i, std: props.permx = try_get_value(*this->fp, "PERMX", props.active_index); props.permy = try_get_value(*this->fp, "PERMY", props.active_index); props.permz = try_get_value(*this->fp, "PERMZ", props.active_index); - //props.satnum = this->fp->int("SATNUM"); + props.satnum = this->fp->get_int("SATNUM").at(props.active_index); + props.pvtnum = this->fp->get_int("PVTNUM").at(props.active_index); + props.ntg = try_get_ntg_value(*this->fp, "NTG", props.active_index); cell.props = props; } } diff --git a/src/opm/parser/eclipse/EclipseState/Schedule/Well/Well.cpp b/src/opm/parser/eclipse/EclipseState/Schedule/Well/Well.cpp index 694eb6732..d5ad70edd 100644 --- a/src/opm/parser/eclipse/EclipseState/Schedule/Well/Well.cpp +++ b/src/opm/parser/eclipse/EclipseState/Schedule/Well/Well.cpp @@ -689,12 +689,12 @@ bool Well::updateConnections(std::shared_ptr connections_arg, b return false; } -bool Well::updateConnections(std::shared_ptr connections_arg, const ScheduleGrid& grid, const std::vector& pvtnum) { +bool Well::updateConnections(std::shared_ptr connections_arg, const ScheduleGrid& grid) { bool update = this->updateConnections(connections_arg, false); if (this->pvt_table == 0 && !this->connections->empty()) { const auto& lowest = this->connections->lowest(); - auto active_index = grid.get_cell(lowest.getI(), lowest.getJ(), lowest.getK()).active_index(); - this->pvt_table = pvtnum[active_index]; + const auto& props = grid.get_cell(lowest.getI(), lowest.getJ(), lowest.getK()).props; + this->pvt_table = props->pvtnum; update = true; } return update; diff --git a/src/opm/parser/eclipse/EclipseState/Schedule/Well/WellConnections.cpp b/src/opm/parser/eclipse/EclipseState/Schedule/Well/WellConnections.cpp index 7a3fba610..4ecdf3bf1 100644 --- a/src/opm/parser/eclipse/EclipseState/Schedule/Well/WellConnections.cpp +++ b/src/opm/parser/eclipse/EclipseState/Schedule/Well/WellConnections.cpp @@ -269,23 +269,8 @@ inline std::array< size_t, 3> directionIndices(const Opm::Connection::Direction defaultSatTabId); } - void WellConnections::loadCOMPDAT(const DeckRecord& record, const ScheduleGrid& grid, const FieldPropsManager& field_properties, const std::string& wname, const KeywordLocation& location) { - const auto permx = field_properties.try_get("PERMX"); - const auto permy = field_properties.try_get("PERMY"); - const auto permz = field_properties.try_get("PERMZ"); - const auto& ntg = field_properties.get_double("NTG"); - const auto& satnum_data = field_properties.get_int("SATNUM"); - - this->loadCOMPDAT(record, grid, satnum_data, permx, permy, permz, ntg, wname, location); - } - void WellConnections::loadCOMPDAT(const DeckRecord& record, const ScheduleGrid& grid, - const std::vector& satnum_data, - const std::vector* permx, - const std::vector* permy, - const std::vector* permz, - const std::vector& ntg, const std::string& wname, const KeywordLocation& location) { @@ -335,15 +320,14 @@ inline std::array< size_t, 3> directionIndices(const Opm::Connection::Direction OpmLog::warning(msg); continue; } - - size_t active_index = cell.active_index(); + const auto& props = cell.props; double CF = -1; double Kh = -1; double r0 = -1; auto ctf_kind = ::Opm::Connection::CTFKind::DeckValue; if (defaultSatTable) - satTableId = satnum_data[active_index]; + satTableId = props->satnum; auto same_ijk = [&]( const Connection& c ) { return c.sameCoordinate( I,J,k ); @@ -362,17 +346,17 @@ inline std::array< size_t, 3> directionIndices(const Opm::Connection::Direction // placement so there's complete exposure (= 2\pi). const double angle = 6.2831853071795864769252867665590057683943387987502116419498; std::array cell_size = cell.dimensions; - const auto& D = effectiveExtent(direction, ntg[active_index], cell_size); + const auto& D = effectiveExtent(direction, props->ntg, cell_size); /* We start with the absolute happy path; both CF and Kh are explicitly given in the deck. */ if (CF > 0 && Kh > 0) goto CF_done; /* We must calculate CF and Kh from the items in the COMPDAT record and cell properties. */ - if (permx && permy && permz) { - std::array cell_perm = {{ permx->operator[](active_index), - permy->operator[](active_index), - permz->operator[](active_index)}}; + { + std::array cell_perm = {{ props->permx, + props->permy, + props->permz}}; const auto& K = permComponents(direction, cell_perm); if (r0 < 0) @@ -391,8 +375,7 @@ inline std::array< size_t, 3> directionIndices(const Opm::Connection::Direction Kh = std::sqrt(K[0] * K[1]) * D[2]; } } - } else - throw std::invalid_argument("Missing PERM values to calculate connection factors"); + } CF_done: if (r0 < 0) diff --git a/tests/msim/compdat.include b/tests/msim/compdat.include new file mode 100644 index 000000000..b4a1d7833 --- /dev/null +++ b/tests/msim/compdat.include @@ -0,0 +1,634 @@ +std::string compdat_deck = R"( +-- This reservoir simulation deck is made available under the Open Database +-- License: http://opendatacommons.org/licenses/odbl/1.0/. Any rights in +-- individual contents of the database are licensed under the Database Contents +-- License: http://opendatacommons.org/licenses/dbcl/1.0/ + +-- Copyright (C) 2015 Statoil + +-- This simulation is based on the data given in +-- 'Comparison of Solutions to a Three-Dimensional +-- Black-Oil Reservoir Simulation Problem' by Aziz S. Odeh, +-- Journal of Petroleum Technology, January 1981 + + +--------------------------------------------------------------------------- +------------------------ SPE1 - CASE 1 ------------------------------------ +--------------------------------------------------------------------------- + +RUNSPEC +-- ------------------------------------------------------------------------- + +TITLE + SPE1 - CASE 1 + +DIMENS + 10 10 3 / + +-- The number of equilibration regions is inferred from the EQLDIMS +-- keyword. +EQLDIMS +/ + +-- The number of PVTW tables is inferred from the TABDIMS keyword; +-- when no data is included in the keyword the default values are used. +TABDIMS +/ + +OIL +GAS +WATER +DISGAS +-- As seen from figure 4 in Odeh, GOR is increasing with time, +-- which means that dissolved gas is present + +FIELD + +START + 1 'DEC' 2014 / + +WELLDIMS +-- Item 1: maximum number of wells in the model +-- - there are two wells in the problem; injector and producer +-- Item 2: maximum number of grid blocks connected to any one well +-- - must be one as the wells are located at specific grid blocks +-- Item 3: maximum number of groups in the model +-- - we are dealing with only one 'group' +-- Item 4: maximum number of wells in any one group +-- - there must be two wells in a group as there are two wells in total + 5 1 1 2 / + +UNIFOUT + +UDQDIMS + 50 25 0 50 50 0 0 50 0 20 / + +GRID + +-- The INIT keyword is used to request an .INIT file. The .INIT file +-- is written before the simulation actually starts, and contains grid +-- properties and saturation tables as inferred from the input +-- deck. There are no other keywords which can be used to configure +-- exactly what is written to the .INIT file. +INIT + + +-- ------------------------------------------------------------------------- +NOECHO + +DX +-- There are in total 300 cells with length 1000ft in x-direction + 300*1000 / +DY +-- There are in total 300 cells with length 1000ft in y-direction + 300*1000 / +DZ +-- The layers are 20, 30 and 50 ft thick, in each layer there are 100 cells + 100*20 100*30 100*50 / + +TOPS +-- The depth of the top of each grid block + 100*8325 / + +PORO +-- Constant porosity of 0.3 throughout all 300 grid cells + 300*0.3 / + +PERMX +-- The layers have perm. 500mD, 50mD and 200mD, respectively. + 100*500 100*50 100*200 / + +PERMY +-- Equal to PERMX + 100*500 100*50 100*200 / + +PERMZ +-- Cannot find perm. in z-direction in Odeh's paper +-- For the time being, we will assume PERMZ equal to PERMX and PERMY: + 100*500 100*50 100*200 / +ECHO + +PROPS +-- ------------------------------------------------------------------------- + +PVTW +-- Item 1: pressure reference (psia) +-- Item 2: water FVF (rb per bbl or rb per stb) +-- Item 3: water compressibility (psi^{-1}) +-- Item 4: water viscosity (cp) +-- Item 5: water 'viscosibility' (psi^{-1}) + +-- Using values from Norne: +-- In METRIC units: +-- 277.0 1.038 4.67E-5 0.318 0.0 / +-- In FIELD units: + 4017.55 1.038 3.22E-6 0.318 0.0 / + +ROCK +-- Item 1: reference pressure (psia) +-- Item 2: rock compressibility (psi^{-1}) + +-- Using values from table 1 in Odeh: + 14.7 3E-6 / + +SWOF +-- Column 1: water saturation +-- - this has been set to (almost) equally spaced values from 0.12 to 1 +-- Column 2: water relative permeability +-- - generated from the Corey-type approx. formula +-- the coeffisient is set to 10e-5, S_{orw}=0 and S_{wi}=0.12 +-- Column 3: oil relative permeability when only oil and water are present +-- - we will use the same values as in column 3 in SGOF. +-- This is not really correct, but since only the first +-- two values are of importance, this does not really matter +-- Column 4: water-oil capillary pressure (psi) + +0.12 0 1 0 +0.18 4.64876033057851E-008 1 0 +0.24 0.000000186 0.997 0 +0.3 4.18388429752066E-007 0.98 0 +0.36 7.43801652892562E-007 0.7 0 +0.42 1.16219008264463E-006 0.35 0 +0.48 1.67355371900826E-006 0.2 0 +0.54 2.27789256198347E-006 0.09 0 +0.6 2.97520661157025E-006 0.021 0 +0.66 3.7654958677686E-006 0.01 0 +0.72 4.64876033057851E-006 0.001 0 +0.78 0.000005625 0.0001 0 +0.84 6.69421487603306E-006 0 0 +0.91 8.05914256198347E-006 0 0 +1 0.00001 0 0 / + + +SGOF +-- Column 1: gas saturation +-- Column 2: gas relative permeability +-- Column 3: oil relative permeability when oil, gas and connate water are present +-- Column 4: oil-gas capillary pressure (psi) +-- - stated to be zero in Odeh's paper + +-- Values in column 1-3 are taken from table 3 in Odeh's paper: +0 0 1 0 +0.001 0 1 0 +0.02 0 0.997 0 +0.05 0.005 0.980 0 +0.12 0.025 0.700 0 +0.2 0.075 0.350 0 +0.25 0.125 0.200 0 +0.3 0.190 0.090 0 +0.4 0.410 0.021 0 +0.45 0.60 0.010 0 +0.5 0.72 0.001 0 +0.6 0.87 0.0001 0 +0.7 0.94 0.000 0 +0.85 0.98 0.000 0 +0.88 0.984 0.000 0 / +--1.00 1.0 0.000 0 / +-- Warning from Eclipse: first sat. value in SWOF + last sat. value in SGOF +-- must not be greater than 1, but Eclipse still runs +-- Flow needs the sum to be excactly 1 so I added a row with gas sat. = 0.88 +-- The corresponding krg value was estimated by assuming linear rel. between +-- gas sat. and krw. between gas sat. 0.85 and 1.00 (the last two values given) + +DENSITY +-- Density (lb per ft³) at surface cond. of +-- oil, water and gas, respectively (in that order) + +-- Using values from Norne: +-- In METRIC units: +-- 859.5 1033.0 0.854 / +-- In FIELD units: + 53.66 64.49 0.0533 / + +PVDG +-- Column 1: gas phase pressure (psia) +-- Column 2: gas formation volume factor (rb per Mscf) +-- - in Odeh's paper the units are said to be given in rb per bbl, +-- but this is assumed to be a mistake: FVF-values in Odeh's paper +-- are given in rb per scf, not rb per bbl. This will be in +-- agreement with conventions +-- Column 3: gas viscosity (cP) + +-- Using values from lower right table in Odeh's table 2: +14.700 166.666 0.008000 +264.70 12.0930 0.009600 +514.70 6.27400 0.011200 +1014.7 3.19700 0.014000 +2014.7 1.61400 0.018900 +2514.7 1.29400 0.020800 +3014.7 1.08000 0.022800 +4014.7 0.81100 0.026800 +5014.7 0.64900 0.030900 +9014.7 0.38600 0.047000 / + +PVTO +-- Column 1: dissolved gas-oil ratio (Mscf per stb) +-- Column 2: bubble point pressure (psia) +-- Column 3: oil FVF for saturated oil (rb per stb) +-- Column 4: oil viscosity for saturated oil (cP) + +-- Use values from top left table in Odeh's table 2: +0.0010 14.7 1.0620 1.0400 / +0.0905 264.7 1.1500 0.9750 / +0.1800 514.7 1.2070 0.9100 / +0.3710 1014.7 1.2950 0.8300 / +0.6360 2014.7 1.4350 0.6950 / +0.7750 2514.7 1.5000 0.6410 / +0.9300 3014.7 1.5650 0.5940 / +1.2700 4014.7 1.6950 0.5100 + 9014.7 1.5790 0.7400 / +1.6180 5014.7 1.8270 0.4490 + 9014.7 1.7370 0.6310 / +-- It is required to enter data for undersaturated oil for the highest GOR +-- (i.e. the last row) in the PVTO table. +-- In order to fulfill this requirement, values for oil FVF and viscosity +-- at 9014.7psia and GOR=1.618 for undersaturated oil have been approximated: +-- It has been assumed that there is a linear relation between the GOR +-- and the FVF when keeping the pressure constant at 9014.7psia. +-- From Odeh we know that (at 9014.7psia) the FVF is 2.357 at GOR=2.984 +-- for saturated oil and that the FVF is 1.579 at GOR=1.27 for undersaturated oil, +-- so it is possible to use the assumption described above. +-- An equivalent approximation for the viscosity has been used. +/ + +SOLUTION +-- ------------------------------------------------------------------------- + +EQUIL +-- Item 1: datum depth (ft) +-- Item 2: pressure at datum depth (psia) +-- - Odeh's table 1 says that initial reservoir pressure is +-- 4800 psi at 8400ft, which explains choice of item 1 and 2 +-- Item 3: depth of water-oil contact (ft) +-- - chosen to be directly under the reservoir +-- Item 4: oil-water capillary pressure at the water oil contact (psi) +-- - given to be 0 in Odeh's paper +-- Item 5: depth of gas-oil contact (ft) +-- - chosen to be directly above the reservoir +-- Item 6: gas-oil capillary pressure at gas-oil contact (psi) +-- - given to be 0 in Odeh's paper +-- Item 7: RSVD-table +-- Item 8: RVVD-table +-- Item 9: Set to 0 as this is the only value supported by OPM + +-- Item #: 1 2 3 4 5 6 7 8 9 + 8400 4800 8450 0 8300 0 1 0 0 / + +RSVD +-- Dissolved GOR is initially constant with depth through the reservoir. +-- The reason is that the initial reservoir pressure given is higher +---than the bubble point presssure of 4014.7psia, meaning that there is no +-- free gas initially present. +8300 1.270 +8450 1.270 / + +SUMMARY +-- ------------------------------------------------------------------------- + +FOPR + +WGOR +/ +WOPR +/ +WWPR +/ +WWCT +/ + + +FGOR + +-- 2a) Pressures of the cell where the injector and producer are located +BPR +1 1 1 / +10 10 3 / +/ + +-- 2b) Gas saturation at grid points given in Odeh's paper +BGSAT +1 1 1 / +1 1 2 / +1 1 3 / +10 1 1 / +10 1 2 / +10 1 3 / +10 10 1 / +10 10 2 / +10 10 3 / +/ + +-- In order to compare Eclipse with Flow: +WBHP +/ + +WGIR + 'INJ' +/ + +WGIT + 'INJ' +/ + +WGPT +/ + +WOPR +/ + +WOPT +/ + +WWIR +/ +WWIT +/ +WWPR +/ +WWPT +/ +WUBHP +/ +WUOPRL +/ +WUWCT +/ + +FOPR + +FUOPR +FUGPR +FU_VAR2 + +SCHEDULE +-- ------------------------------------------------------------------------- +RPTSCHED + 'PRES' 'SGAS' 'RS' 'WELLS' 'WELSPECS' / + +RPTRST + 'BASIC=1' / + +UDQ + ASSIGN WUBHP 11 / + ASSIGN WUOPRL 20 / + ASSIGN WUBHP P2 12 / + ASSIGN WUBHP P3 13 / + ASSIGN WUBHP P4 14 / + UNITS WUBHP 'BARSA' / + UNITS WUOPRL 'SM3/DAY' / + DEFINE WUWCT WWPR / (WWPR + WOPR) / + UNITS WUWCT '1' / + DEFINE FUOPR SUM(WOPR) / + UNITS FUOPR 'SM3/DAY' / + DEFINE FUGPR FLPR / +/ + + + +-- If no resolution (i.e. case 1), the two following lines must be added: +DRSDT + 0 / +-- if DRSDT is set to 0, GOR cannot rise and free gas does not +-- dissolve in undersaturated oil -> constant bubble point pressure + +WELSPECS +-- Item #: 1 2 3 4 5 6 + 'P1' 'G1' 3 3 8400 'OIL' / + 'P2' 'G1' 4 4 8400 'OIL' / + 'P3' 'G1' 5 5 8400 'OIL' / + 'P4' 'G1' 6 6 8400 'OIL' / + 'INJ' 'G1' 1 1 8335 'GAS' / +/ +-- Coordinates in item 3-4 are retrieved from Odeh's figure 1 and 2 +-- Note that the depth at the midpoint of the well grid blocks +-- has been used as reference depth for bottom hole pressure in item 5 + +COMPDAT +-- Item #: 1 2 3 4 5 6 7 8 9 + 'P1' 3 3 3 3 'OPEN' 1* 1* 0.5 / + 'P2' 4 4 3 3 'OPEN' 1* 1* 0.5 / + 'P3' 5 5 3 3 'OPEN' 1* 1* 0.5 / + 'P4' 6 6 3 3 'OPEN' 1* 1* 0.5 / + 'INJ' 1 1 1 1 'OPEN' 1* 1* 0.5 / +/ +-- Coordinates in item 2-5 are retreived from Odeh's figure 1 and 2 +-- Item 9 is the well bore internal diameter, +-- the radius is given to be 0.25ft in Odeh's paper + +UDQ + DEFINE FUINJ 0.90 * SUM(WWPR 'P*') / + UNITS FUINJ 'SM3/DAY' / +/ + + + +ACTIONX + 'SHUT_WELL' 100000 / + WWCT * > 0.50 / +/ + +UDQ + DEFINE FUNEW 0.90 * SUM(WWPR 'P*') / + UNITS FUNEW 'SM3/DAY' / +/ + +WELOPEN + '?' 'SHUT' 0 0 0 2* / +/ + +ENDACTIO + + +WCONPROD +-- Item #:1 2 3 4 5 9 + 'P1' 'OPEN' 'ORAT' 5000 4* 1000 / + 'P2' 'OPEN' 'ORAT' 5000 4* 1000 / + 'P3' 'OPEN' 'ORAT' 5000 4* 1000 / + 'P4' 'OPEN' 'ORAT' 5000 4* 1000 / +/ + +-- It is stated in Odeh's paper that the maximum oil prod. rate +-- is 20 000stb per day which explains the choice of value in item 4. +-- The items > 4 are defaulted with the exception of item 9, +-- the BHP lower limit, which is given to be 1000psia in Odeh's paper + +WCONINJE +-- Item #:1 2 3 4 5 6 7 + 'INJ' 'GAS' 'OPEN' 'RATE' 100000 1* 9014 / +/ +-- Stated in Odeh that gas inj. rate (item 5) is 100MMscf per day +-- BHP upper limit (item 7) should not be exceeding the highest +-- pressure in the PVT table=9014.7psia (default is 100 000psia) + +DATES -- 1 + 1 'JAN' 2015 / +/ + +DATES -- 2 + 1 'FEB' 2015 / +/ + +DATES -- 3 + 1 'MAR' 2015 / +/ + +UDQ + DEFINE FU_VAR2 FGLIR / +/ + + + +DATES -- 4 + 1 'APR' 2015 / +/ + +DATES -- 5 + 1 'MAI' 2015 / +/ + +UDQ + DEFINE FU_VAR3 WGPR P1 / +/ + +DATES + 1 'JUN' 2015 / +/ + +DATES + 1 'JUL' 2015 / +/ + +DATES + 1 'AUG' 2015 / +/ + +DATES + 1 'SEP' 2015 / +/ + +DATES -- 10 + 1 'OCT' 2015 / +/ + +DATES + 1 'NOV' 2015 / +/ + +DATES + 1 'DEC' 2015 / +/ + +DATES + 1 'JAN' 2016 / +/ + +DATES + 1 'FEB' 2016 / +/ + +DATES -- 15 + 1 'MAR' 2016 / +/ + +DATES + 1 'APR' 2016 / +/ + +DATES + 1 'MAI' 2016 / +/ + +DATES + 1 'JUN' 2016 / +/ + +DATES + 1 'JUL' 2016 / +/ + +DATES -- 20 + 1 'AUG' 2016 / +/ + +DATES + 1 'SEP' 2016 / +/ + +DATES + 1 'OCT' 2016 / +/ + +DATES + 1 'NOV' 2016 / +/ + +DATES -- 24 + 1 'DEC' 2016 / +/ + +UDQ + DEFINE FUPROD 0.90 * SUM(WWPR 'P*') / + UNITS FUPROD 'SM3/DAY' / +/ + + + +DATES -- 25 + 1 'JAN' 2017 / +/ + +DATES + 1 'FEB' 2017 / +/ + +DATES + 1 'MAR' 2017 / +/ + +DATES + 1 'APR' 2017 / +/ + +DATES + 1 'MAI' 2017 / +/ + +DATES -- 30 + 1 'JUN' 2017 / +/ + +DATES + 1 'JUL' 2017 / +/ + +DATES + 1 'AUG' 2017 / +/ + +DATES + 1 'SEP' 2017 / +/ + +DATES + 1 'OCT' 2017 / +/ + +DATES -- 35 + 1 'NOV' 2017 / +/ + +COMPDAT +-- Item #: 1 2 3 4 5 6 7 8 9 + 'P4' 6 6 1 3 'OPEN' 1* 1* 0.5 / +/ + + + +DATES -- 36 + 1 'DEC' 2017 / +/ + + + +END +)"; diff --git a/tests/msim/test_msim_ACTIONX.cpp b/tests/msim/test_msim_ACTIONX.cpp index 5fe8d22bb..cfe328c8e 100644 --- a/tests/msim/test_msim_ACTIONX.cpp +++ b/tests/msim/test_msim_ACTIONX.cpp @@ -457,6 +457,24 @@ BOOST_AUTO_TEST_CASE(UDA) { } } +BOOST_AUTO_TEST_CASE(COMPDAT) { +#include "compdat.include" + test_data td( compdat_deck ); + msim sim(td.state); + EclipseIO io(td.state, td.state.getInputGrid(), td.schedule, td.summary_config); + + sim.well_rate("P1", data::Rates::opt::wat, prod_wpr_P1); + sim.well_rate("P2", data::Rates::opt::wat, prod_wpr_P2); + sim.well_rate("P3", data::Rates::opt::wat, prod_wpr_P3); + sim.well_rate("P4", data::Rates::opt::wat, prod_wpr_P4); + sim.well_rate("INJ", data::Rates::opt::wat, inj_wir_INJ); + { + WorkArea work_area("compdat_sim"); + + BOOST_CHECK_NO_THROW(sim.run(td.schedule, io, true)); + } +} + #ifdef EMBEDDED_PYTHON BOOST_AUTO_TEST_CASE(PYTHON_WELL_CLOSE_EXAMPLE) { diff --git a/tests/parser/ACTIONX.cpp b/tests/parser/ACTIONX.cpp index 8d278a20b..5667101af 100644 --- a/tests/parser/ACTIONX.cpp +++ b/tests/parser/ACTIONX.cpp @@ -185,6 +185,49 @@ TSTEP BOOST_CHECK_THROW( make_schedule(WITH_GRID, parseContext), OpmInputError ); } +BOOST_AUTO_TEST_CASE(COMPDAT) { + + const auto TRAILING_COMPDAT = std::string{ R"( +GRID + +PORO + 1000*0.1 / +PERMX + 1000*1 / +PERMY + 1000*0.1 / +PERMZ + 1000*0.01 / + +SCHEDULE + +WELSPECS + 'W2' 'OP' 1 1 3.33 'OIL' 7*/ +/ + +ACTIONX + 'ACTION' / + WWCT OPX > 0.75 / +/ + +ENDACTIO + +TSTEP + 10 / + +COMPDAT + 'W2' 1 1 1 1 'OPEN' / +/ + +)"}; + + Schedule sched = make_schedule(TRAILING_COMPDAT); + Action::Result action_result(true); + auto sim_time = TimeService::now(); + const auto& action1 = sched[0].actions.get()["ACTION"]; + BOOST_CHECK_NO_THROW( sched.applyAction(0, sim_time, action1, Action::Result{true}, {})); +} + BOOST_AUTO_TEST_CASE(TestActions) { Opm::SummaryState st(TimeService::now()); diff --git a/tests/parser/ConnectionTests.cpp b/tests/parser/ConnectionTests.cpp index 0d9b30049..201a2cd88 100644 --- a/tests/parser/ConnectionTests.cpp +++ b/tests/parser/ConnectionTests.cpp @@ -59,7 +59,7 @@ Opm::WellConnections loadCOMPDAT(const std::string& compdat_keyword) { Opm::WellConnections connections(Opm::Connection::Order::TRACK, 10,10); Opm::CompletedCells cells(grid); for (const auto& rec : keyword) - connections.loadCOMPDAT(rec, Opm::ScheduleGrid(grid, field_props, cells), field_props, "WELL", {}); + connections.loadCOMPDAT(rec, Opm::ScheduleGrid(grid, field_props, cells), "WELL", {}); return connections; }