OilfieldToolsNet/opm-simulators
4
0
Fork 0
mirror of https://github.com/OPM/opm-simulators.git synced 2026-09-05 04:40:19 -05:00
Code Issues Packages Projects Releases Wiki Activity

Merge pull request #4180 from bska/report-segment-dyn-data

Browse Source 
Convey More Dynamic State Values per Segment to Output Layer
This commit is contained in:
Bård Skaflestad
2022-11-17 12:39:29 +01:00
committed by GitHub
parent 3c31b7bafe 936e9e87e7
commit 4845ad0107
7 changed files with 263 additions and 75 deletions
Download Patch File Download Diff File Expand all files Collapse all files
opm/simulators/wells/MultisegmentWellEval.cpp
+125 -8
View File
@@ -19,6 +19,7 @@
*/
#include <config.h>
#include <opm/simulators/wells/MultisegmentWellEval.hpp>
#include <dune/istl/umfpack.hh>
@@ -33,13 +34,20 @@
#include <opm/simulators/timestepping/ConvergenceReport.hpp>
#include <opm/simulators/utils/DeferredLoggingErrorHelpers.hpp>
#include <opm/simulators/wells/MSWellHelpers.hpp>
#include <opm/simulators/wells/RateConverter.hpp>
#include <opm/simulators/wells/WellAssemble.hpp>
#include <opm/simulators/wells/WellBhpThpCalculator.hpp>
#include <opm/simulators/wells/WellConvergence.hpp>
#include <opm/simulators/wells/WellInterfaceIndices.hpp>
#include <opm/simulators/wells/WellState.hpp>
#include <algorithm>
#include <array>
#include <cassert>
#include <cmath>
#include <numeric>
#include <utility>
#include <vector>
namespace Opm
{
@@ -1397,6 +1405,8 @@ updateWellStateFromPrimaryVariables(WellState& well_state,
static constexpr int Oil = BlackoilPhases::Liquid;
static constexpr int Water = BlackoilPhases::Aqua;
const auto pvtReg = std::max(this->baseif_.wellEcl().pvt_table_number() - 1, 0);
const PhaseUsage& pu = baseif_.phaseUsage();
assert( FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) );
const int oil_pos = pu.phase_pos[Oil];
@@ -1409,13 +1419,13 @@ updateWellStateFromPrimaryVariables(WellState& well_state,
std::vector<double> fractions(baseif_.numPhases(), 0.0);
fractions[oil_pos] = 1.0;
if ( FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx) ) {
if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
const int water_pos = pu.phase_pos[Water];
fractions[water_pos] = primary_variables_[seg][WFrac];
fractions[oil_pos] -= fractions[water_pos];
}
if ( FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx) ) {
if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
const int gas_pos = pu.phase_pos[Gas];
fractions[gas_pos] = primary_variables_[seg][GFrac];
fractions[oil_pos] -= fractions[gas_pos];
@@ -1445,16 +1455,123 @@ updateWellStateFromPrimaryVariables(WellState& well_state,
// update the segment pressure
segment_pressure[seg] = primary_variables_[seg][SPres];
if (seg == 0) { // top segment
ws.bhp = segment_pressure[seg];
}
// Calculate other per-phase dynamic quantities.
const auto temperature = 0.0; // Ignore thermal effects
const auto saltConc = 0.0; // Ignore salt precipitation
const auto Rvw = 0.0; // Ignore vaporised water.
auto rsMax = 0.0;
auto rvMax = 0.0;
if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
// Both oil and gas active.
rsMax = FluidSystem::oilPvt()
.saturatedGasDissolutionFactor(pvtReg, temperature, segment_pressure[seg]);
rvMax = FluidSystem::gasPvt()
.saturatedOilVaporizationFactor(pvtReg, temperature, segment_pressure[seg]);
}
// 1) Local condition volume flow rates
const auto& [Rs, Rv] = this->baseif_.rateConverter().inferDissolvedVaporisedRatio
(rsMax, rvMax, segment_rates.begin() + (seg + 0)*this->baseif_.numPhases());
{
// Use std::span<> in C++20 and beyond.
const auto rate_start = (seg + 0) * this->baseif_.numPhases();
const auto* surf_rates = segment_rates.data() + rate_start;
auto* resv_rates = segments.phase_resv_rates.data() + rate_start;
this->baseif_.rateConverter().calcReservoirVoidageRates
(pvtReg, segment_pressure[seg],
std::max(0.0, Rs),
std::max(0.0, Rv),
temperature, saltConc, surf_rates, resv_rates);
}
// 2) Local condition holdup fractions.
const auto tot_resv =
std::accumulate(segments.phase_resv_rates.begin() + (seg + 0)*this->baseif_.numPhases(),
segments.phase_resv_rates.begin() + (seg + 1)*this->baseif_.numPhases(),
0.0);
std::transform(segments.phase_resv_rates.begin() + (seg + 0)*this->baseif_.numPhases(),
segments.phase_resv_rates.begin() + (seg + 1)*this->baseif_.numPhases(),
segments.phase_holdup.begin() + (seg + 0)*this->baseif_.numPhases(),
[tot_resv](const auto qr) { return std::clamp(qr / tot_resv, 0.0, 1.0); });
// 3) Local condition flow velocities for segments other than top segment.
if (seg > 0) {
// Possibly poor approximation
// Velocity = Flow rate / cross-sectional area.
// Additionally ignores drift flux.
const auto area = this->baseif_.wellEcl().getSegments()
.getFromSegmentNumber(segments.segment_number()[seg]).crossArea();
const auto velocity = (area > 0.0) ? tot_resv / area : 0.0;
std::transform(segments.phase_holdup.begin() + (seg + 0)*this->baseif_.numPhases(),
segments.phase_holdup.begin() + (seg + 1)*this->baseif_.numPhases(),
segments.phase_velocity.begin() + (seg + 0)*this->baseif_.numPhases(),
[velocity](const auto hf) { return (hf > 0.0) ? velocity : 0.0; });
}
// 4) Local condition phase viscosities.
segments.phase_viscosity[seg*this->baseif_.numPhases() + pu.phase_pos[Oil]] =
FluidSystem::oilPvt().viscosity(pvtReg, temperature, segment_pressure[seg], Rs);
if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
segments.phase_viscosity[seg*this->baseif_.numPhases() + pu.phase_pos[Water]] =
FluidSystem::waterPvt().viscosity(pvtReg, temperature, segment_pressure[seg], saltConc);
}
if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
segments.phase_viscosity[seg*this->baseif_.numPhases() + pu.phase_pos[Gas]] =
FluidSystem::gasPvt().viscosity(pvtReg, temperature, segment_pressure[seg], Rv, Rvw);
}
}
WellBhpThpCalculator(baseif_).
updateThp(rho, [this]() { return baseif_.wellEcl().alq_value(); },
{FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx),
FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx),
FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)},
well_state, deferred_logger);
// Segment flow velocity in top segment.
{
const auto np = this->baseif_.numPhases();
auto segVel = [&segments, np](const auto segmentNumber)
{
auto v = 0.0;
const auto* vel = segments.phase_velocity.data() + segmentNumber*np;
for (auto p = 0*np; p < np; ++p) {
if (std::abs(vel[p]) > std::abs(v)) {
v = vel[p];
}
}
return v;
};
const auto seg = 0;
auto maxVel = 0.0;
for (const auto& inlet : this->segmentSet()[seg].inletSegments()) {
const auto v = segVel(this->segmentNumberToIndex(inlet));
if (std::abs(v) > std::abs(maxVel)) {
maxVel = v;
}
}
std::transform(segments.phase_holdup.begin() + (seg + 0)*this->baseif_.numPhases(),
segments.phase_holdup.begin() + (seg + 1)*this->baseif_.numPhases(),
segments.phase_velocity.begin() + (seg + 0)*this->baseif_.numPhases(),
[maxVel](const auto hf) { return (hf > 0.0) ? maxVel : 0.0; });
}
WellBhpThpCalculator(this->baseif_)
.updateThp(rho, [this]() { return this->baseif_.wellEcl().alq_value(); },
{FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx),
FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx),
FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)},
well_state, deferred_logger);
}
template<typename FluidSystem, typename Indices, typename Scalar>
opm/simulators/wells/MultisegmentWellEval.hpp
+2
View File
@@ -35,6 +35,8 @@
#include <array>
#include <memory>
#include <utility>
#include <vector>
namespace Dune {
template<class Matrix> class UMFPack;
opm/simulators/wells/RateConverter.hpp
+22 -11
View File
@@ -422,18 +422,16 @@ namespace Opm {
* \param[out] voidage_rates Reservoir volume flow rates for all
* active phases.
*/
template <class Rates>
void calcReservoirVoidageRates(const int pvtRegionIdx,
const double p,
const double rs,
const double rv,
const double T,
const double saltConcentration,
const Rates& surface_rates,
Rates& voidage_rates) const
template <typename SurfaceRates, typename VoidageRates>
void calcReservoirVoidageRates(const int pvtRegionIdx,
const double p,
const double rs,
const double rv,
const double T,
const double saltConcentration,
const SurfaceRates& surface_rates,
VoidageRates& voidage_rates) const
{
std::fill(voidage_rates.begin(), voidage_rates.end(), 0.0);
const auto& pu = this->phaseUsage_;
const auto iw = RegionAttributeHelpers::PhasePos::water(pu);
const auto io = RegionAttributeHelpers::PhasePos::oil (pu);
@@ -442,6 +440,8 @@ namespace Opm {
const auto [Rs, Rv] = this->
dissolvedVaporisedRatio(io, ig, rs, rv, surface_rates);
std::fill_n(&voidage_rates[0], pu.num_phases, 0.0);
if (RegionAttributeHelpers::PhaseUsed::water(pu)) {
// q[w]_r = q[w]_s / bw
const auto bw = FluidSystem::waterPvt()
@@ -482,6 +482,17 @@ namespace Opm {
}
}
template <class Rates>
std::pair<double, double>
inferDissolvedVaporisedRatio(const double rsMax,
const double rvMax,
const Rates& surface_rates) const
{
const auto io = RegionAttributeHelpers::PhasePos::oil(this->phaseUsage_);
const auto ig = RegionAttributeHelpers::PhasePos::gas(this->phaseUsage_);
return this->dissolvedVaporisedRatio(io, ig, rsMax, rvMax, surface_rates);
}
/**
* Compute coefficients for surface-to-reservoir voidage
* conversion for solvent.
opm/simulators/wells/SegmentState.cpp
+39 -25
View File
@@ -16,45 +16,60 @@
You should have received a copy of the GNU General Public License
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#if HAVE_CONFIG_H
#include "config.h"
#endif // HAVE_CONFIG_H
#include <algorithm>
#include <stdexcept>
#include <opm/simulators/wells/SegmentState.hpp>
#include <opm/input/eclipse/Schedule/Well/WellConnections.hpp>
#include <opm/input/eclipse/Schedule/MSW/WellSegments.hpp>
#include <algorithm>
#include <cstddef>
#include <iterator>
#include <stdexcept>
#include <vector>
namespace {
std::vector<int> make_segment_number(const Opm::WellSegments& segments)
{
std::vector<int> segment_number;
segment_number.reserve(segments.size());
std::transform(segments.begin(), segments.end(),
std::back_inserter(segment_number),
[](const Opm::Segment& segment)
{
return segment.segmentNumber();
});
return segment_number;
}
} // Anonymous namespace
namespace Opm
{
namespace {
std::vector<int> make_segment_number( const WellSegments& segments ) {
std::vector<int> segment_number;
std::transform(segments.begin(), segments.end(), std::back_insert_iterator(segment_number), [](const Segment& segment) { return segment.segmentNumber(); });
return segment_number;
}
}
SegmentState::SegmentState(int num_phases, const WellSegments& segments) :
rates(segments.size() * num_phases),
pressure(segments.size()),
pressure_drop_friction(segments.size()),
pressure_drop_hydrostatic(segments.size()),
pressure_drop_accel(segments.size()),
m_segment_number(make_segment_number(segments))
{
}
SegmentState::SegmentState(int num_phases, const WellSegments& segments)
: rates (segments.size() * num_phases)
, phase_resv_rates (segments.size() * num_phases)
, phase_velocity (segments.size() * num_phases)
, phase_holdup (segments.size() * num_phases)
, phase_viscosity (segments.size() * num_phases)
, pressure (segments.size())
, pressure_drop_friction (segments.size())
, pressure_drop_hydrostatic(segments.size())
, pressure_drop_accel (segments.size())
, m_segment_number (make_segment_number(segments))
{}
double SegmentState::pressure_drop(std::size_t index) const {
return this->pressure_drop_friction[index] + this->pressure_drop_hydrostatic[index] + this->pressure_drop_accel[index];
}
bool SegmentState::empty() const {
return this->rates.empty();
}
@@ -63,7 +78,6 @@ std::size_t SegmentState::size() const {
return this->pressure.size();
}
void SegmentState::scale_pressure(const double bhp) {
if (this->empty())
throw std::logic_error("Tried to pressure scale empty SegmentState");
@@ -80,4 +94,4 @@ const std::vector<int>& SegmentState::segment_number() const {
return this->m_segment_number;
}
}
} // namespace Opm
opm/simulators/wells/SegmentState.hpp
+22 -4
View File
@@ -20,35 +20,53 @@
#ifndef OPM_SEGMENTSTATE_HEADER_INCLUDED
#define OPM_SEGMENTSTATE_HEADER_INCLUDED
#include <cstddef>
#include <vector>
namespace Opm
{
class WellSegments;
class WellConnections;
} // namespace Opm
namespace Opm
{
class SegmentState
{
public:
SegmentState() = default;
SegmentState(int num_phases, const WellSegments& segments);
double pressure_drop(std::size_t index) const;
bool empty() const;
void scale_pressure(double bhp);
const std::vector<int>& segment_number() const;
std::size_t size() const;
std::vector<double> rates;
/// Segment condition volume flow rates through segment (per phase)
std::vector<double> phase_resv_rates;
/// Segment condition flow velocity through segment (per phase)
std::vector<double> phase_velocity;
/// Segment condition holdup fractions through segment (per phase)
std::vector<double> phase_holdup;
/// Segment condition phase viscosities.
std::vector<double> phase_viscosity;
std::vector<double> pressure;
std::vector<double> pressure_drop_friction;
std::vector<double> pressure_drop_hydrostatic;
std::vector<double> pressure_drop_accel;
private:
std::vector<int> m_segment_number;
};
}
} // namepace Opm
#endif
#endif // OPM_SEGMENTSTATE_HEADER_INCLUDED
opm/simulators/wells/WellState.cpp
+45 -21
View File
@@ -18,15 +18,15 @@
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#include <fmt/format.h>
#include <config.h>
#include <opm/simulators/wells/WellState.hpp>
#include <opm/common/ErrorMacros.hpp>
#include <opm/input/eclipse/Schedule/Schedule.hpp>
#include <opm/simulators/wells/ParallelWellInfo.hpp>
#include <opm/simulators/utils/ParallelCommunication.hpp>
#include <opm/simulators/wells/ParallelWellInfo.hpp>
#include <opm/grid/common/p2pcommunicator.hh>
#include <opm/output/data/Wells.hpp>
@@ -35,6 +35,9 @@
#include <numeric>
#include <set>
#include <stdexcept>
#include <vector>
#include <fmt/format.h>
namespace {
@@ -425,14 +428,15 @@ data::Wells
WellState::report(const int* globalCellIdxMap,
const std::function<bool(const int)>& wasDynamicallyClosed) const
{
if (this->numWells() == 0)
if (this->numWells() == 0) {
return {};
}
using rt = data::Rates::opt;
const auto& pu = this->phaseUsage();
data::Wells res;
for( std::size_t well_index = 0; well_index < this->size(); well_index++) {
for (std::size_t well_index = 0; well_index < this->size(); ++well_index) {
const auto& ws = this->well(well_index);
if ((ws.status == Well::Status::SHUT) && !wasDynamicallyClosed(well_index))
{
@@ -505,12 +509,10 @@ WellState::report(const int* globalCellIdxMap,
}
const auto& pwinfo = ws.parallel_info.get();
if (pwinfo.communication().size()==1)
{
if (pwinfo.communication().size() == 1) {
reportConnections(well.connections, pu, well_index, globalCellIdxMap);
}
else
{
else {
std::vector<data::Connection> connections;
reportConnections(connections, pu, well_index, globalCellIdxMap);
@@ -828,17 +830,21 @@ void WellState::updateGlobalIsGrup(const Comm& comm)
}
data::Segment
WellState::reportSegmentResults(const int well_id,
const int seg_ix,
const int seg_no) const
WellState::reportSegmentResults(const int well_id,
const int seg_ix,
const int seg_no) const
{
using PhaseQuant = data::SegmentPhaseQuantity::Item;
const auto& segments = this->well(well_id).segments;
if (segments.empty())
if (segments.empty()) {
return {};
}
auto seg_res = data::Segment{};
{
using Value = data::SegmentPressures::Value;
auto& segpress = seg_res.pressures;
segpress[Value::Pressure] = segments.pressure[seg_ix];
segpress[Value::PDrop] = segments.pressure_drop(seg_ix);
@@ -848,20 +854,40 @@ WellState::reportSegmentResults(const int well_id,
}
const auto& pu = this->phaseUsage();
const auto rate = &segments.rates[seg_ix * pu.num_phases];
const auto* rate = &segments.rates[seg_ix * pu.num_phases];
const auto* resv = &segments.phase_resv_rates[seg_ix * pu.num_phases];
const auto* velocity = &segments.phase_velocity[seg_ix * pu.num_phases];
const auto* holdup = &segments.phase_holdup[seg_ix * pu.num_phases];
const auto* viscosity = &segments.phase_viscosity[seg_ix * pu.num_phases];
if (pu.phase_used[Water]) {
seg_res.rates.set(data::Rates::opt::wat,
rate[pu.phase_pos[Water]]);
const auto iw = pu.phase_pos[Water];
seg_res.rates.set(data::Rates::opt::wat, rate[iw]);
seg_res.rates.set(data::Rates::opt::reservoir_water, resv[iw]);
seg_res.velocity.set(PhaseQuant::Water, velocity[iw]);
seg_res.holdup.set(PhaseQuant::Water, holdup[iw]);
seg_res.viscosity.set(PhaseQuant::Water, viscosity[iw]);
}
if (pu.phase_used[Oil]) {
seg_res.rates.set(data::Rates::opt::oil,
rate[pu.phase_pos[Oil]]);
const auto io = pu.phase_pos[Oil];
seg_res.rates.set(data::Rates::opt::oil, rate[io]);
seg_res.rates.set(data::Rates::opt::reservoir_oil, resv[io]);
seg_res.velocity.set(PhaseQuant::Oil, velocity[io]);
seg_res.holdup.set(PhaseQuant::Oil, holdup[io]);
seg_res.viscosity.set(PhaseQuant::Oil, viscosity[io]);
}
if (pu.phase_used[Gas]) {
seg_res.rates.set(data::Rates::opt::gas,
rate[pu.phase_pos[Gas]]);
const auto ig = pu.phase_pos[Gas];
seg_res.rates.set(data::Rates::opt::gas, rate[ig]);
seg_res.rates.set(data::Rates::opt::reservoir_gas, resv[ig]);
seg_res.velocity.set(PhaseQuant::Gas, velocity[ig]);
seg_res.holdup.set(PhaseQuant::Gas, holdup[ig]);
seg_res.viscosity.set(PhaseQuant::Gas, viscosity[ig]);
}
seg_res.segNumber = seg_no;
@@ -913,5 +939,3 @@ WellState::parallelWellInfo(std::size_t well_index) const
template void WellState::updateGlobalIsGrup<Parallel::Communication>(const Parallel::Communication& comm);
template void WellState::communicateGroupRates<Parallel::Communication>(const Parallel::Communication& comm);
} // namespace Opm
opm/simulators/wells/WellState.hpp
+8 -6
View File
@@ -21,14 +21,16 @@
#ifndef OPM_WELLSTATEFULLYIMPLICITBLACKOIL_HEADER_INCLUDED
#define OPM_WELLSTATEFULLYIMPLICITBLACKOIL_HEADER_INCLUDED
#include <opm/simulators/wells/ALQState.hpp>
#include <opm/simulators/wells/SingleWellState.hpp>
#include <opm/simulators/wells/GlobalWellInfo.hpp>
#include <opm/simulators/wells/SegmentState.hpp>
#include <opm/simulators/wells/WellContainer.hpp>
#include <opm/core/props/BlackoilPhases.hpp>
#include <opm/simulators/wells/PerforationData.hpp>
#include <opm/simulators/wells/ALQState.hpp>
#include <opm/simulators/wells/GlobalWellInfo.hpp>
#include <opm/simulators/wells/PerfData.hpp>
#include <opm/simulators/wells/PerforationData.hpp>
#include <opm/simulators/wells/SegmentState.hpp>
#include <opm/simulators/wells/SingleWellState.hpp>
#include <opm/simulators/wells/WellContainer.hpp>
#include <opm/output/data/Wells.hpp>
#include <opm/input/eclipse/Schedule/Events.hpp>