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

Merge pull request #1576 from joakim-hove/rst-msw-init

Browse Source 
Rst msw init
This commit is contained in:
Joakim Hove 2020-03-13 15:45:18 +01:00 committed by GitHub
commit 7e88f4b4c2
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
6 changed files with 63 additions and 11 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

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

@ -32,6 +32,7 @@ class Header;
struct RstConnection {
RstConnection(const ::Opm::UnitSystem& unit_system, 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::array<int,3> ijk;

4
opm/parser/eclipse/EclipseState/Schedule/TimeMap.hpp
View File

@ -27,6 +27,7 @@
#include <ctime>
#include <map>
#include <utility>
#include <iostream>
#include <stddef.h>
@ -99,6 +100,9 @@ namespace Opm {
bool m_skiprest = false;
std::size_t m_restart_offset = 0;
};
std::ostream& operator<<(std::ostream& stream, const TimeMap& tm);
}

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

@ -66,7 +66,11 @@ Connection::CTFKind from_float(float float_kind) {
return Connection::CTFKind::DeckValue;
}
}
double RstConnection::inverse_peaceman(double cf, double kh, double rw, double skin) {
auto alpha = 3.14159265 * 2 * kh / cf - skin;
return rw * std::exp(alpha);
}
using M = ::Opm::UnitSystem::measure;
@ -92,11 +96,19 @@ RstConnection::RstConnection(const ::Opm::UnitSystem& unit_system, const int* ic
water_rate( unit_system.to_si(M::liquid_surface_rate, xcon[VI::XConn::WaterRate])),
gas_rate( unit_system.to_si(M::gas_surface_rate, xcon[VI::XConn::GasRate])),
pressure( unit_system.to_si(M::pressure, xcon[VI::XConn::Pressure])),
resv_rate( unit_system.to_si(M::rate, xcon[VI::XConn::ResVRate]))
{
auto alpha = 3.14159265 * 2 * this->kh / this->cf - this->skin_factor;
this->r0 = this->diameter * std::exp(alpha) / 2;
}
resv_rate( unit_system.to_si(M::rate, xcon[VI::XConn::ResVRate])),
r0( RstConnection::inverse_peaceman(this->cf, this->kh, this->diameter/2, this->skin_factor) )
/*
r0: The r0 quantity is currently not written or read from the restart
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.
*/
{ }
}
}

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

@ -3120,6 +3120,19 @@ bool Schedule::cmp(const Schedule& sched1, const Schedule& sched2, std::size_t r
if (count != 0)
return false;
{
const auto& tm1 = sched1.getTimeMap();
const auto& tm2 = sched2.getTimeMap();
if (not_equal(tm1.size(), tm2.size(), "TimeMap: size()"))
count += 1;
for (auto& step_index = report_step; step_index < std::min(tm1.size(), tm2.size()) - 1; step_index++) {
if (not_equal(tm1[step_index], tm2[step_index], "TimePoint[" + std::to_string(step_index) + "]"))
count += 1;
}
}
for (const auto& wname : sched1.wellNames(report_step)) {
const auto& well1 = sched1.getWell(wname, report_step);
const auto& well2 = sched2.getWell(wname, report_step);

18
src/opm/parser/eclipse/EclipseState/Schedule/TimeMap.cpp
View File

@ -20,6 +20,7 @@
#include <cassert>
#include <ctime>
#include <stddef.h>
#include <iomanip>
#include <opm/common/utility/TimeService.hpp>
@ -436,6 +437,23 @@ namespace {
bool TimeMap::skiprest() const {
return this->m_skiprest;
}
std::ostream& operator<<(std::ostream& stream, const TimeMap& tm) {
std::stringstream ss;
ss << "{";
std::size_t index = 0;
for (const auto& tp : tm.timeList()) {
auto ts = TimeStampUTC(tp);
ss << ts.year() << "-" << std::setfill('0') << std::setw(2) << ts.month() << "-" << std::setfill('0') << std::setw(2) << ts.day();
index += 1;
if (index < tm.timeList().size())
ss << ", ";
if (index % 12 == 0)
ss << std::endl;
}
ss << "}";
return stream << ss.str();
}
}

16
src/opm/parser/eclipse/EclipseState/Schedule/Well/WellConnections.cpp
View File

@ -23,6 +23,7 @@
#include <iostream>
#include <opm/parser/eclipse/Units/Units.hpp>
#include <opm/io/eclipse/rst/connection.hpp>
#include <opm/parser/eclipse/EclipseState/Grid/EclipseGrid.hpp>
#include <opm/parser/eclipse/EclipseState/Grid/FieldPropsManager.hpp>
#include <opm/parser/eclipse/EclipseState/Schedule/Well/Connection.hpp>
@ -244,16 +245,14 @@ inline std::array< size_t, 3> directionIndices(const Opm::Connection::Direction
int satTableId = -1;
bool defaultSatTable = true;
const auto& r0Item = record.getItem("PR");
const auto& CFItem = record.getItem("CONNECTION_TRANSMISSIBILITY_FACTOR");
const auto& diameterItem = record.getItem("DIAMETER");
const auto& KhItem = record.getItem("Kh");
const auto& satTableIdItem = record.getItem("SAT_TABLE");
const auto& r0Item = record.getItem("PR");
const auto direction = Connection::DirectionFromString(record.getItem("DIR").getTrimmedString(0));
double skin_factor = record.getItem("SKIN").getSIDouble(0);
double rw;
double r0=0.0;
if (satTableIdItem.hasValue(0) && satTableIdItem.get < int > (0) > 0)
{
@ -276,6 +275,7 @@ inline std::array< size_t, 3> directionIndices(const Opm::Connection::Direction
size_t active_index = grid.activeIndex(I,J,k);
double CF = -1;
double Kh = -1;
double r0 = -1;
auto ctf_kind = ::Opm::Connection::CTFKind::DeckValue;
if (defaultSatTable)
@ -285,6 +285,9 @@ inline std::array< size_t, 3> directionIndices(const Opm::Connection::Direction
return c.sameCoordinate( I,J,k );
};
if (r0Item.hasValue(0))
r0 = r0Item.getSIDouble(0);
if (KhItem.hasValue(0) && KhItem.getSIDouble(0) > 0.0)
Kh = KhItem.getSIDouble(0);
@ -307,9 +310,7 @@ inline std::array< size_t, 3> directionIndices(const Opm::Connection::Direction
const auto& K = permComponents(direction, cell_perm);
const auto& D = effectiveExtent(direction, ntg[active_index], cell_size);
if (r0Item.hasValue(0))
r0 = r0Item.getSIDouble(0);
else
if (r0 < 0)
r0 = effectiveRadius(K,D);
if (CF < 0) {
@ -330,6 +331,9 @@ inline std::array< size_t, 3> directionIndices(const Opm::Connection::Direction
CF_done:
if (r0 < 0)
r0 = RestartIO::RstConnection::inverse_peaceman(CF, Kh, rw, skin_factor);
auto prev = std::find_if( this->m_connections.begin(),
this->m_connections.end(),
same_ijk );