Add functions for well rates and solutions to mock simulator

This commit is contained in:
Joakim Hove 2018-11-15 14:31:41 +01:00
parent 6af3194c6c
commit bdcb227a53
3 changed files with 108 additions and 3 deletions

View File

@ -1,6 +1,10 @@
#ifndef ISIM_MAIN_HPP
#define ISIM_MAIN_HPP
#include <functional>
#include <string>
#include <map>
#include <opm/parser/eclipse/Deck/Deck.hpp>
#include <opm/parser/eclipse/EclipseState/EclipseState.hpp>
#include <opm/parser/eclipse/EclipseState/EclipseState.hpp>
@ -15,22 +19,33 @@ namespace Opm {
class EclipseIO;
class msim {
public:
using well_rate_function = double(const EclipseState&, const Schedule&, const data::Solution&, size_t report_step, double seconds_elapsed);
using solution_function = void(const EclipseState&, const Schedule&, data::Solution&, size_t report_step, double seconds_elapsed);
msim(const std::string& deck_file);
msim(const std::string& deck_file, const Parser& parser, const ParseContext& parse_context);
void well_rate(const std::string& well, data::Rates::opt rate, std::function<well_rate_function> func);
void solution(const std::string& field, std::function<solution_function> func);
void run();
private:
void run_step(data::Solution& sol, data::Wells& well_data, size_t report_step, EclipseIO& io) const;
void run_step(data::Solution& sol, data::Wells& well_data, size_t report_step, double dt, EclipseIO& io) const;
void output(size_t report_step, bool substep, double seconds_elapsed, const data::Solution& sol, const data::Wells& well_data, EclipseIO& io) const;
void simulate(data::Solution& sol, data::Wells& well_data, size_t report_step, double seconds_elapsed, double time_step) const;
Deck deck;
EclipseState state;
Schedule schedule;
SummaryConfig summary_config;
std::map<std::string, std::map<data::Rates::opt, std::function<well_rate_function>>> well_rates;
std::map<std::string, std::function<solution_function>> solutions;
};
}

View File

@ -20,7 +20,7 @@
#include <iostream>
#include <opm/output/eclipse/EclipseIO.hpp>
#include <opm/output/eclipse/RestartValue.hpp>
#include <opm/output/eclipse/RestartValue.hpp>
#include <opm/output/data/Solution.hpp>
#include <opm/output/data/Wells.hpp>
@ -73,7 +73,7 @@ void msim::run_step(data::Solution& sol, data::Wells& well_data, size_t report_s
if ((seconds_elapsed + time_step) > end_time)
time_step = end_time - seconds_elapsed;
// Simulate
this->simulate(sol, well_data, report_step, seconds_elapsed, time_step);
seconds_elapsed += time_step;
this->output(report_step,
@ -98,4 +98,36 @@ void msim::output(size_t report_step, bool substep, double seconds_elapsed, cons
{});
}
void msim::simulate(data::Solution& sol, data::Wells& well_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, this->schedule, sol, report_step, seconds_elapsed + time_step);
}
for (const auto& well_pair : this->well_rates) {
const std::string& well_name = well_pair.first;
data::Well& well = well_data[well_name];
for (const auto& rate_pair : well_pair.second) {
auto rate = rate_pair.first;
auto func = rate_pair.second;
well.rates.set(rate, func(this->state, this->schedule, sol, report_step, seconds_elapsed + time_step));
}
}
}
void msim::well_rate(const std::string& well, data::Rates::opt rate, std::function<well_rate_function> func) {
this->well_rates[well][rate] = func;
}
void msim::solution(const std::string& field, std::function<solution_function> func) {
this->solutions[field] = func;
}
}

View File

@ -26,21 +26,79 @@
#include <ert/util/test_work_area.h>
#include <ert/util/util.h>
#include <ert/ecl/ecl_sum.hpp>
#include <ert/ecl/ecl_file.hpp>
#include <ert/ecl/ecl_file_view.hpp>
#include <ert/ecl/ecl_rsthead.hpp>
#include <opm/msim/msim.hpp>
#include <opm/output/data/Wells.hpp>
#include <opm/parser/eclipse/Units/Units.hpp>
using namespace Opm;
double prod_opr(const EclipseState& es, const Schedule& sched, const data::Solution& sol, size_t report_step, double seconds_elapsed) {
const auto& units = es.getUnits();
return -units.to_si(UnitSystem::measure::rate, seconds_elapsed);
}
void pressure(const EclipseState& es, const Schedule& sched, data::Solution& sol, size_t report_step, double seconds_elapsed) {
const auto& grid = es.getInputGrid();
const auto& units = es.getUnits();
if (!sol.has("PRESSURE"))
sol.insert("PRESSURE", UnitSystem::measure::pressure, std::vector<double>(grid.getNumActive()), data::TargetType::RESTART_SOLUTION);
auto& data = sol.data("PRESSURE");
std::fill(data.begin(), data.end(), units.to_si(UnitSystem::measure::pressure, seconds_elapsed));
}
BOOST_AUTO_TEST_CASE(RUN) {
msim msim("SPE1CASE1.DATA");
msim.well_rate("PROD", data::Rates::opt::oil, prod_opr);
msim.solution("PRESSURE", pressure);
{
test_work_area_type * work_area = test_work_area_alloc("test_msim");
msim.run();
for (const auto& fname : {"SPE1CASE1.INIT", "SPE1CASE1.UNRST", "SPE1CASE1.EGRID", "SPE1CASE1.SMSPEC"})
for (const auto& fname : {"SPE1CASE1.INIT", "SPE1CASE1.UNRST", "SPE1CASE1.EGRID", "SPE1CASE1.SMSPEC", "SPE1CASE1.UNSMRY"})
BOOST_CHECK( util_is_file( fname ));
{
ecl_sum_type * ecl_sum = ecl_sum_fread_alloc_case("SPE1CASE1", ":");
int param_index = ecl_sum_get_general_var_params_index(ecl_sum, "WOPR:PROD");
for (int time_index=0; time_index < ecl_sum_get_data_length(ecl_sum); time_index++) {
double seconds_elapsed = ecl_sum_iget_sim_days(ecl_sum, time_index) * 86400;
double opr = ecl_sum_iget(ecl_sum, time_index, param_index);
opr = ecl_sum_iget(ecl_sum, time_index, param_index);
BOOST_CHECK_CLOSE(seconds_elapsed, opr, 1e-3);
}
ecl_sum_free( ecl_sum );
}
{
int report_step = 1;
ecl_file_type * rst_file = ecl_file_open("SPE1CASE1.UNRST", 0);
ecl_file_view_type * global_view = ecl_file_get_global_view( rst_file );
while (true) {
ecl_file_view_type * rst_view = ecl_file_view_add_restart_view( global_view, -1, report_step, -1, -1);
if (!rst_view)
break;
{
ecl_rsthead_type * rst_head = ecl_rsthead_alloc( rst_view, report_step);
const ecl_kw_type * p = ecl_file_view_iget_named_kw(rst_view, "PRESSURE", 0);
BOOST_CHECK_CLOSE( ecl_kw_iget_float(p, 0), rst_head->sim_days * 86400, 1e-3 );
ecl_rsthead_free( rst_head );
}
report_step++;
}
ecl_file_close(rst_file);
}
test_work_area_free( work_area );
}
}