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

Merge pull request #3173 from bska/write-rft-seg-data

Browse Source 
Enable SEG Type RFT File Output
This commit is contained in:
Markus Blatt 2022-11-22 12:23:32 +01:00 committed by GitHub
commit 5590b3c400
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
4 changed files with 1731 additions and 73 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

1
opm/input/eclipse/Schedule/MSW/WellSegments.hpp
View File

@ -74,6 +74,7 @@ namespace Opm {
static WellSegments serializationTestObject();
std::size_t size() const;
bool empty() const;
double depthTopSegment() const;
double lengthTopSegment() const;
double volumeTopSegment() const;

4
src/opm/input/eclipse/Schedule/MSW/WellSegments.cpp
View File

@ -73,6 +73,10 @@ namespace Opm {
return m_segments.size();
}
bool WellSegments::empty() const {
return this->m_segments.empty();
}
const Segment& WellSegments::topSegment() const {
return this->m_segments[0];
}

542
src/opm/output/eclipse/WriteRFT.cpp
View File

@ -53,6 +53,7 @@
#include <optional>
#include <stdexcept>
#include <string>
#include <type_traits>
#include <utility>
#include <vector>
@ -364,7 +365,60 @@ namespace {
// =======================================================================
class PLTFlowRate
class PLTPhaseQuantity
{
public:
explicit PLTPhaseQuantity(const std::size_t n = 0);
virtual ~PLTPhaseQuantity() = default;
const std::vector<float>& oil() const { return this->oil_; }
const std::vector<float>& gas() const { return this->gas_; }
const std::vector<float>& water() const { return this->water_; }
protected:
void addOil (const ::Opm::UnitSystem& usys, const float x);
void addGas (const ::Opm::UnitSystem& usys, const float x);
void addWater(const ::Opm::UnitSystem& usys, const float x);
private:
std::vector<float> oil_{};
std::vector<float> gas_{};
std::vector<float> water_{};
[[nodiscard]] virtual ::Opm::UnitSystem::measure oilUnit() const = 0;
[[nodiscard]] virtual ::Opm::UnitSystem::measure gasUnit() const = 0;
[[nodiscard]] virtual ::Opm::UnitSystem::measure waterUnit() const = 0;
};
void PLTPhaseQuantity::addOil(const ::Opm::UnitSystem& usys, const float x)
{
this->oil_.push_back(usys.from_si(this->oilUnit(), x));
}
void PLTPhaseQuantity::addGas(const ::Opm::UnitSystem& usys, const float x)
{
this->gas_.push_back(usys.from_si(this->gasUnit(), x));
}
void PLTPhaseQuantity::addWater(const ::Opm::UnitSystem& usys, const float x)
{
this->water_.push_back(usys.from_si(this->waterUnit(), x));
}
PLTPhaseQuantity::PLTPhaseQuantity(const std::size_t n)
{
if (n == std::size_t{0}) {
return;
}
this->oil_.reserve(n);
this->gas_.reserve(n);
this->water_.reserve(n);
}
// -----------------------------------------------------------------------
class PLTFlowRate : public PLTPhaseQuantity
{
public:
explicit PLTFlowRate(const std::size_t nconn = 0);
@ -372,40 +426,33 @@ namespace {
void addConnection(const ::Opm::UnitSystem& usys,
const ::Opm::data::Rates& rates);
const std::vector<float>& oil() const { return this->oil_; }
const std::vector<float>& gas() const { return this->gas_; }
const std::vector<float>& water() const { return this->water_; }
private:
std::vector<float> oil_{};
std::vector<float> gas_{};
std::vector<float> water_{};
[[nodiscard]] Opm::UnitSystem::measure oilUnit() const override
{ return Opm::UnitSystem::measure::liquid_surface_rate; }
[[nodiscard]] Opm::UnitSystem::measure gasUnit() const override
{ return Opm::UnitSystem::measure::gas_surface_rate; }
[[nodiscard]] Opm::UnitSystem::measure waterUnit() const override
{ return Opm::UnitSystem::measure::liquid_surface_rate; }
};
PLTFlowRate::PLTFlowRate(const std::size_t nconn)
{
if (nconn == std::size_t{0}) {
return;
}
this->oil_.reserve(nconn);
this->gas_.reserve(nconn);
this->water_.reserve(nconn);
}
: PLTPhaseQuantity{nconn}
{}
void PLTFlowRate::addConnection(const ::Opm::UnitSystem& usys,
const ::Opm::data::Rates& rates)
{
using M = Opm::UnitSystem::measure;
using rt = Opm::data::Rates::opt;
// Note negative sign on call to rates.get() here. Flow reports
// positive injection rates and negative production rates but we
// need the opposite sign convention for this report.
this->oil_.push_back(usys.from_si(M::liquid_surface_rate, -rates.get(rt::oil, 0.0)));
this->gas_.push_back(usys.from_si(M::gas_surface_rate, -rates.get(rt::gas, 0.0)));
this->water_.push_back(usys.from_si(M::liquid_surface_rate, -rates.get(rt::wat, 0.0)));
this->addOil (usys, -rates.get(rt::oil, 0.0));
this->addGas (usys, -rates.get(rt::gas, 0.0));
this->addWater(usys, -rates.get(rt::wat, 0.0));
}
// -----------------------------------------------------------------------
@ -566,6 +613,8 @@ namespace {
const std::function<int(int)>& binId,
Cmp&& cmp);
int maxBinId() const { return this->maxId_; }
IndexRange bin(const int binId) const
{
this->verifyValid(binId);
@ -679,6 +728,54 @@ namespace {
// -----------------------------------------------------------------------
class OrderSegments
{
public:
explicit OrderSegments(const ::Opm::WellSegments& wellSegs)
: wellSegs_{ std::cref(wellSegs) }
{}
bool operator()(const int i1, const int i2) const;
private:
std::reference_wrapper<const ::Opm::WellSegments> wellSegs_;
};
// i1 < i2 if one of the following relations hold
//
// 1) i1's branch number is smaller than i2's branch number
// 2) i1 and i2 are on the same branch, but i1 is i2's outlet segment
// 3) Neither are each other's outlet segments, but i1 is closer to the
// well head along the tubing.
bool OrderSegments::operator()(const int i1, const int i2) const
{
const auto& s1 = this->wellSegs_.get()[i1];
const auto& s2 = this->wellSegs_.get()[i2];
const auto b1 = s1.branchNumber();
const auto b2 = s2.branchNumber();
if (b1 != b2) {
// i1 not on same branch as i2. Order by branch number.
return b1 < b2;
}
if (s2.outletSegment() == s1.segmentNumber()) {
// i1 is i2's outlet
return true;
}
if (s1.outletSegment() == s2.segmentNumber()) {
// i2 is i1's outlet
return false;
}
// Neither is each other's outlet. Order by distance along tubing.
return s1.totalLength() < s2.totalLength();
}
// -----------------------------------------------------------------------
class PLTRecordMSW : public PLTRecord
{
public:
@ -687,19 +784,6 @@ namespace {
void write(::Opm::EclIO::OutputStream::RFT& rftFile) const override;
private:
class OrderSegments
{
public:
explicit OrderSegments(const ::Opm::WellSegments& wellSegs)
: wellSegs_{ std::cref(wellSegs) }
{}
bool operator()(const int i1, const int i2) const;
private:
std::reference_wrapper<const ::Opm::WellSegments> wellSegs_;
};
class OrderSegConns
{
public:
@ -880,41 +964,6 @@ namespace {
// -----------------------------------------------------------------------
// i1 < i2 if one of the following relations hold
//
// 1) i1's branch number is smaller than i2's branch number
// 2) i1 and i2 are on the same branch, but i1 is i2's outlet segment
// 3) Neither are each other's outlet segments, but i1 is closer to the
// well head along the tubing.
bool PLTRecordMSW::OrderSegments::operator()(const int i1, const int i2) const
{
const auto& s1 = this->wellSegs_.get()[i1];
const auto& s2 = this->wellSegs_.get()[i2];
const auto b1 = s1.branchNumber();
const auto b2 = s2.branchNumber();
if (b1 != b2) {
// i1 not on same branch as i2. Order by branch number.
return b1 < b2;
}
if (s2.outletSegment() == s1.segmentNumber()) {
// i1 is i2's outlet
return true;
}
if (s1.outletSegment() == s2.segmentNumber()) {
// i2 is i1's outlet
return false;
}
// Neither is each other's outlet. Order by distance along tubing.
return s1.totalLength() < s2.totalLength();
}
// -----------------------------------------------------------------------
// i1 < i2 if one of the following relations hold
//
// 1) i1's branch number is smaller than i2's branch number
@ -964,13 +1013,320 @@ namespace {
return this->wellConns_.get()[connIdx].perf_range()->second;
}
// -----------------------------------------------------------------------
class SegmentRecord
{
public:
explicit SegmentRecord(const std::size_t nseg = 0);
void collectRecordData(const ::Opm::UnitSystem& usys,
const ::Opm::Well& well,
const ::Opm::data::Well& wellSol);
std::size_t nSeg() const { return this->neighbour_id_.size(); }
void write(::Opm::EclIO::OutputStream::RFT& rftFile) const;
private:
PLTFlowRate rate_{};
std::vector<int> neighbour_id_{};
std::vector<int> branch_id_{};
std::vector<int> branch_start_segment_{};
std::vector<int> branch_end_segment_{};
std::vector<float> diameter_{};
std::vector<float> depth_{};
std::vector<float> start_length_{};
std::vector<float> end_length_{};
std::vector<float> x_coord_{};
std::vector<float> y_coord_{};
std::vector<float> pressure_{};
std::vector<float> strength_{};
std::vector<float> icd_setting_{};
void defineBranches(const ::Opm::WellSegments& segments);
void addSegment(const ::Opm::UnitSystem& usys,
const ::Opm::WellSegments& segments,
const ::Opm::Segment& segment,
const ::Opm::data::Segment& segSol);
void recordPhysicalLocation(const ::Opm::UnitSystem& usys,
const ::Opm::WellSegments& segments,
const ::Opm::Segment& segment);
void recordSegmentConnectivity(const ::Opm::Segment& segment);
void recordSegmentProperties(const ::Opm::UnitSystem& usys,
const ::Opm::Segment& segment);
void recordDynamicState(const ::Opm::UnitSystem& usys,
const ::Opm::Segment& segment,
const ::Opm::data::Segment& segSol);
void recordAutoICDTypeProperties(const ::Opm::UnitSystem& usys,
const ::Opm::Segment& segment);
void recordRegularTypeProperties(const ::Opm::UnitSystem& usys,
const ::Opm::Segment& segment);
void recordSpiralICDTypeProperties(const ::Opm::UnitSystem& usys,
const ::Opm::Segment& segment);
void recordValveTypeProperties(const ::Opm::UnitSystem& usys,
const ::Opm::Segment& segment);
};
SegmentRecord::SegmentRecord(const std::size_t nseg)
: rate_{ nseg }
{
if (nseg == std::size_t{0}) {
return;
}
this->neighbour_id_.reserve(nseg);
this->branch_id_.reserve(nseg);
this->diameter_.reserve(nseg);
this->depth_.reserve(nseg);
this->start_length_.reserve(nseg);
this->end_length_.reserve(nseg);
this->x_coord_.reserve(nseg);
this->y_coord_.reserve(nseg);
this->pressure_.reserve(nseg);
this->strength_.reserve(nseg);
this->icd_setting_.reserve(nseg);
}
void SegmentRecord::write(::Opm::EclIO::OutputStream::RFT& rftFile) const
{
rftFile.write("SEGDIAM", this->diameter_);
rftFile.write("SEGDEPTH", this->depth_);
rftFile.write("SEGLENST", this->start_length_);
rftFile.write("SEGLENEN", this->end_length_);
rftFile.write("SEGXCORD", this->x_coord_);
rftFile.write("SEGYCORD", this->y_coord_);
rftFile.write("SEGPRES", this->pressure_);
rftFile.write("SEGORAT", this->rate_.oil());
rftFile.write("SEGWRAT", this->rate_.water());
rftFile.write("SEGGRAT", this->rate_.gas());
rftFile.write("SEGSSTR", this->strength_);
rftFile.write("SEGSFOPN", this->icd_setting_);
rftFile.write("SEGBRNO", this->branch_id_);
rftFile.write("SEGNXT", this->neighbour_id_);
rftFile.write("BRNST", this->branch_start_segment_);
rftFile.write("BRNEN", this->branch_end_segment_);
}
void SegmentRecord::collectRecordData(const ::Opm::UnitSystem& usys,
const ::Opm::Well& well,
const ::Opm::data::Well& wellSol)
{
const auto& segments = well.getSegments();
const auto& xseg = wellSol.segments;
for (const auto& segment : segments) {
auto segSolPos = xseg.find(segment.segmentNumber());
if (segSolPos == xseg.end()) {
continue;
}
this->addSegment(usys, segments, segment, segSolPos->second);
}
if (this->nSeg() > std::size_t{0}) {
this->defineBranches(segments);
}
}
void SegmentRecord::defineBranches(const ::Opm::WellSegments& wellSegs)
{
auto branchSegments = CSRIndexRelation {};
const auto minBranchID = 1;
branchSegments.build(wellSegs.size(), minBranchID,
[&wellSegs](const int ix) { return wellSegs[ix].branchNumber(); },
OrderSegments { wellSegs });
const auto maxBranchID = branchSegments.maxBinId();
this->branch_start_segment_.assign(maxBranchID, 0);
this->branch_end_segment_ .assign(maxBranchID, 0);
auto segNum = [&wellSegs](const int segIx) {
return wellSegs[segIx].segmentNumber();
};
for (auto branch = minBranchID; branch <= maxBranchID; ++branch) {
auto [begin, end] = branchSegments.bin(branch);
if (end == begin) {
// Empty branch (no segments). Unexpected.
continue;
}
const auto branchIx = branch - minBranchID;
this->branch_start_segment_[branchIx] = segNum(*begin);
this->branch_end_segment_ [branchIx] = segNum(*(end - 1));
}
}
void SegmentRecord::addSegment(const ::Opm::UnitSystem& usys,
const ::Opm::WellSegments& segments,
const ::Opm::Segment& segment,
const Opm::data::Segment& segSol)
{
this->recordPhysicalLocation(usys, segments, segment);
this->recordSegmentConnectivity(segment);
this->recordSegmentProperties(usys, segment);
this->recordDynamicState(usys, segment, segSol);
}
void SegmentRecord::recordPhysicalLocation(const ::Opm::UnitSystem& usys,
const ::Opm::WellSegments& segments,
const ::Opm::Segment& segment)
{
auto cvrt = [&usys](const double x) -> float
{
return usys.from_si(::Opm::UnitSystem::measure::length, x);
};
{
const auto diam = std::max(0.0, segment.internalDiameter());
this->diameter_.push_back(cvrt(diam));
}
this->depth_.push_back(cvrt(segment.depth()));
{
const auto outletNum = segment.outletSegment();
const auto start = (outletNum <= 0)
? 0.0 // 'segment' is top segment
: segments.getFromSegmentNumber(outletNum).totalLength();
this->start_length_.push_back(cvrt(start));
}
this->end_length_.push_back(cvrt(segment.totalLength()));
this->x_coord_.push_back(cvrt(segment.node_X()));
this->y_coord_.push_back(cvrt(segment.node_Y()));
}
void SegmentRecord::recordSegmentConnectivity(const ::Opm::Segment& segment)
{
this->neighbour_id_.push_back(segment.outletSegment());
this->branch_id_ .push_back(segment.branchNumber());
}
void SegmentRecord::recordSegmentProperties(const ::Opm::UnitSystem& usys,
const ::Opm::Segment& segment)
{
using TypeProperties = void(SegmentRecord::*)
(const ::Opm::UnitSystem& usys,
const ::Opm::Segment& segment);
static const auto handlers = std::array<TypeProperties, 4> {
&SegmentRecord::recordRegularTypeProperties,
&SegmentRecord::recordSpiralICDTypeProperties,
&SegmentRecord::recordAutoICDTypeProperties,
&SegmentRecord::recordValveTypeProperties,
};
const auto ix = static_cast<std::size_t>(segment.segmentType());
if (ix < handlers.size()) {
(this->*handlers[ix])(usys, segment);
}
else {
// Unknown segment type. Treat as "regular".
this->recordRegularTypeProperties(usys, segment);
}
}
void SegmentRecord::recordDynamicState(const ::Opm::UnitSystem& usys,
[[maybe_unused]] const ::Opm::Segment& segment,
const ::Opm::data::Segment& segSol)
{
using M = ::Opm::UnitSystem::measure;
using SegPress = ::Opm::data::SegmentPressures::Value;
this->pressure_.push_back(usys.from_si(M::pressure, segSol.pressures[SegPress::Pressure]));
this->rate_.addConnection(usys, segSol.rates);
}
void SegmentRecord::recordAutoICDTypeProperties(const ::Opm::UnitSystem& usys,
const ::Opm::Segment& segment)
{
this->strength_.push_back(usys.from_si(::Opm::UnitSystem::measure::aicd_strength,
segment.autoICD().strength()));
this->icd_setting_.push_back(1.0f);
}
void SegmentRecord::recordRegularTypeProperties([[maybe_unused]] const ::Opm::UnitSystem& usys,
[[maybe_unused]] const ::Opm::Segment& segment)
{
this->strength_.push_back(0.0f);
this->icd_setting_.push_back(1.0f);
}
void SegmentRecord::recordSpiralICDTypeProperties(const ::Opm::UnitSystem& usys,
const ::Opm::Segment& segment)
{
this->strength_.push_back(usys.from_si(::Opm::UnitSystem::measure::icd_strength,
segment.spiralICD().strength()));
this->icd_setting_.push_back(1.0f);
}
double valveMaximumCrossSectionalArea(const ::Opm::Segment& segment)
{
assert ((segment.segmentNumber() > 1) && segment.isValve());
// Data sources in order of preference:
//
// 1. WSEGVALV(10) (= valve.conMaxCrossArea()), if set
// 2. WSEGVALV( 8) (= valve.pipeCrossArea()), if set
// 3. WELSEGS ( 9) (= segment.crossArea())
const auto& valve = segment.valve();
if (const auto ac_max = valve.conMaxCrossArea(); ac_max > 0.0) {
return ac_max;
}
if (const auto ac_max = valve.pipeCrossArea(); ac_max > 0.0) {
return ac_max;
}
return segment.crossArea();
}
// ICD setting ("SEGSFOPN") for valves is cross-sectional area in valve
// constriction (Ac) relative to maximum cross-sectional area in value
// constriction (Ac_max).
double valveICDSetting(const ::Opm::Segment& segment)
{
assert ((segment.segmentNumber() > 1) && segment.isValve());
const auto Ac = segment.valve().conCrossArea();
const auto Ac_max = valveMaximumCrossSectionalArea(segment);
return Ac / Ac_max;
}
void SegmentRecord::recordValveTypeProperties([[maybe_unused]] const ::Opm::UnitSystem& usys,
const ::Opm::Segment& segment)
{
this->strength_.push_back(0.0f);
this->icd_setting_.push_back(valveICDSetting(segment));
}
// =======================================================================
class WellRFTOutputData
{
public:
enum class DataTypes {
RFT, PLT,
RFT, PLT, SEG,
};
explicit WellRFTOutputData(const std::vector<DataTypes>& types,
@ -1008,6 +1364,7 @@ namespace {
std::unique_ptr<WellConnectionRecord> wconns_{};
std::unique_ptr<RFTRecord> rft_{};
std::unique_ptr<PLTRecord> plt_{};
std::unique_ptr<SegmentRecord> seg_{};
std::vector<DataHandler> dataHandlers_{};
std::vector<RecordWriter> recordWriters_{};
@ -1017,10 +1374,12 @@ namespace {
void initialiseConnHandlers();
void initialiseRFTHandlers();
void initialisePLTHandlers();
void initialiseSEGHandlers();
bool haveOutputData() const;
bool haveRFTData() const;
bool havePLTData() const;
bool haveSEGData() const;
void writeHeader(::Opm::EclIO::OutputStream::RFT& rftFile) const;
@ -1056,7 +1415,8 @@ namespace {
bool WellRFTOutputData::haveOutputData() const
{
return this->haveRFTData()
|| this->havePLTData();
|| this->havePLTData()
|| this->haveSEGData();
}
void WellRFTOutputData::addDynamicData(const Opm::data::Well& wellSol)
@ -1150,6 +1510,30 @@ namespace {
});
}
void WellRFTOutputData::initialiseSEGHandlers()
{
const auto& well = this->well_.get();
if (!well.isMultiSegment() || well.getSegments().empty()) {
return;
}
this->seg_ = std::make_unique<SegmentRecord>
(well.getSegments().size());
this->dataHandlers_.emplace_back(
[this](const Opm::data::Well& wellSol)
{
this->seg_->collectRecordData(this->usys_, this->well_, wellSol);
});
this->recordWriters_.emplace_back(
[this](::Opm::EclIO::OutputStream::RFT& rftFile)
{
this->seg_->write(rftFile);
});
}
bool WellRFTOutputData::haveRFTData() const
{
return (this->rft_ != nullptr)
@ -1162,6 +1546,12 @@ namespace {
&& (this->plt_->nConn() > std::size_t{0});
}
bool WellRFTOutputData::haveSEGData() const
{
return (this->seg_ != nullptr)
&& (this->seg_->nSeg() > std::size_t{0});
}
void WellRFTOutputData::writeHeader(::Opm::EclIO::OutputStream::RFT& rftFile) const
{
{
@ -1229,8 +1619,9 @@ namespace {
std::string WellRFTOutputData::dataTypeString() const
{
auto tstring = std::string{};
if (this->haveRFTData()) { tstring += 'R';}
if (this->havePLTData()) { tstring += 'P';}
if (this->haveRFTData()) { tstring += 'R'; }
if (this->havePLTData()) { tstring += 'P'; }
if (this->haveSEGData()) { tstring += 'S'; }
return tstring;
}
@ -1246,6 +1637,7 @@ namespace {
WellRFTOutputData::creators_{
{ WellRFTOutputData::DataTypes::RFT, &WellRFTOutputData::initialiseRFTHandlers },
{ WellRFTOutputData::DataTypes::PLT, &WellRFTOutputData::initialisePLTHandlers },
{ WellRFTOutputData::DataTypes::SEG, &WellRFTOutputData::initialiseSEGHandlers },
};
} // Anonymous namespace
@ -1266,6 +1658,10 @@ namespace {
rftTypes.push_back(WellRFTOutputData::DataTypes::PLT);
}
if (rft_config.segment(well_name)) {
rftTypes.push_back(WellRFTOutputData::DataTypes::SEG);
}
return rftTypes;
}
}

1257
tests/test_RFT.cpp
View File

File diff suppressed because it is too large Load Diff