Merge pull request #2908 from bska/dont-access-missing-refdepth

Don't Access Well BHP Reference Depth If Unavailable
This commit is contained in:
Joakim Hove 2021-12-19 12:13:29 +01:00 committed by GitHub
commit d867d63b78
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
6 changed files with 190 additions and 73 deletions

View File

@ -526,6 +526,7 @@ public:
int getHeadI() const;
int getHeadJ() const;
double getWPaveRefDepth() const;
bool hasRefDepth() const;
double getRefDepth() const;
double getDrainageRadius() const;
double getEfficiencyFactor() const;

View File

@ -27,6 +27,26 @@
#include <opm/output/eclipse/VectorItems/well.hpp>
#include <opm/parser/eclipse/Units/UnitSystem.hpp>
#include <cmath>
namespace {
bool is_sentinel(const float raw_value)
{
const auto infty = 1.0e+20f;
return ! (std::abs(raw_value) < infty);
}
double swel_value(const float raw_value)
{
return is_sentinel(raw_value) ? 0.0 : raw_value;
}
template <typename Convert>
double keep_sentinel(const float raw_value, Convert&& convert)
{
return is_sentinel(raw_value) ? raw_value : convert(raw_value);
}
}
namespace VI = ::Opm::RestartIO::Helpers::VectorItems;
@ -38,14 +58,6 @@ constexpr int def_pvt_table = 0;
using M = ::Opm::UnitSystem::measure;
double swel_value(float raw_value) {
const auto infty = 1.0e+20f;
if (std::abs(raw_value) == infty)
return 0;
else
return raw_value;
}
RstWell::RstWell(const ::Opm::UnitSystem& unit_system,
const RstHeader& header,
const std::string& group_arg,
@ -89,7 +101,7 @@ RstWell::RstWell(const ::Opm::UnitSystem& unit_system,
hist_lrat_target( unit_system.to_si(M::liquid_surface_rate, swel[VI::SWell::HistLiqRateTarget])),
hist_grat_target( unit_system.to_si(M::gas_surface_rate, swel[VI::SWell::HistGasRateTarget])),
hist_bhp_target( unit_system.to_si(M::pressure, swel[VI::SWell::HistBHPTarget])),
datum_depth( unit_system.to_si(M::length, swel[VI::SWell::DatumDepth])),
datum_depth( keep_sentinel(swel[VI::SWell::DatumDepth], [&unit_system](const double depth) { return unit_system.to_si(M::length, depth); })),
drainage_radius( unit_system.to_si(M::length, swel_value(swel[VI::SWell::DrainageRadius]))),
efficiency_factor( unit_system.to_si(M::identity, swel[VI::SWell::EfficiencyFactor1])),
alq_value( swel[VI::SWell::Alq_value]),

View File

@ -54,6 +54,7 @@
#include <exception>
#include <iterator>
#include <iostream>
#include <optional>
#include <stdexcept>
#include <string>
@ -435,7 +436,7 @@ namespace {
return inteHead[VI::intehead::NSWELZ];
}
float datumDepth(const Opm::Well& well)
std::optional<float> datumDepth(const Opm::Well& well)
{
if (well.isMultiSegment()) {
// Datum depth for multi-segment wells is
@ -443,9 +444,14 @@ namespace {
return well.getSegments().depthTopSegment();
}
// Not a multi-segment well--i.e., this is a regular
// well. Use well's reference depth.
return well.getRefDepth();
// Not a multi-segment well--i.e., this is a regular well. Use
// well's reference depth if available. Otherwise signal a
// missing value by std::nullopt. Missing values *might* come
// from WELSPECS(5) being defaulted in wells with no active
// reservoir connections.
return well.hasRefDepth()
? std::optional<float>{ well.getRefDepth() }
: std::nullopt;
}
Opm::RestartIO::Helpers::WindowedArray<float>
@ -508,7 +514,9 @@ namespace {
std::copy(b, e, std::begin(sWell));
}
float getRateLimit(const Opm::UnitSystem& units, Opm::UnitSystem::measure u, const double& rate)
float getRateLimit(const Opm::UnitSystem& units,
const Opm::UnitSystem::measure u,
const double rate)
{
float rLimit = 1.0e+20f;
if (rate > 0.0) {
@ -521,7 +529,6 @@ namespace {
return rLimit;
}
template <class SWellArray>
void assignTracerData(const Opm::TracerConfig& tracers,
const Opm::SummaryState& smry,
@ -539,14 +546,14 @@ namespace {
}
template <class SWellArray>
void staticContrib(const Opm::Well& well,
const Opm::GasLiftOpt& glo,
const std::size_t sim_step,
const Opm::Schedule& sched,
const Opm::TracerConfig& tracers,
const Opm::WellTestState& wtest_state,
const ::Opm::SummaryState& smry,
SWellArray& sWell)
void staticContrib(const Opm::Well& well,
const Opm::GasLiftOpt& glo,
const std::size_t sim_step,
const Opm::Schedule& sched,
const Opm::TracerConfig& tracers,
const Opm::WellTestState& wtest_state,
const ::Opm::SummaryState& smry,
SWellArray& sWell)
{
using Ix = ::Opm::RestartIO::Helpers::VectorItems::SWell::index;
using M = ::Opm::UnitSystem::measure;
@ -694,15 +701,27 @@ namespace {
sWell[Ix::LOincFac] = static_cast<float>(w_glo.inc_weight_factor());
}
sWell[Ix::DatumDepth] = swprop(M::length, datumDepth(well));
if (const auto depth = datumDepth(well); depth.has_value()) {
sWell[Ix::DatumDepth] = swprop(M::length, depth.value());
}
else {
// BHP reference depth missing for this well. Typically
// caused by defaulted WELSPECS(5) and no active reservoir
// connections from which to infer the depth. Output
// sentinel value.
//
// Note: All unit systems get the *same* sentinel value so
// we intentionally omit unit conversion here.
sWell[Ix::DatumDepth] = -1.0e+20f;
}
sWell[Ix::DrainageRadius] = swprop(M::length, well.getDrainageRadius());
/*
Restart files from Eclipse indicate that the efficiency factor is
found in two items in the restart file; since only one of the
items is needed in the OPM restart, and we are not really certain
that the two values are equal by construction we have only
assigned one of the items explicitly here.
*/
// Restart files from Eclipse indicate that the efficiency
// factor is found in two items in the restart file; since only
// one of the items is needed in the OPM restart, and we are not
// really certain that the two values are equal by construction
// we have only assigned one of the items explicitly here.
sWell[Ix::EfficiencyFactor1] = well.getEfficiencyFactor();
sWell[Ix::EfficiencyFactor2] = sWell[Ix::EfficiencyFactor1];

View File

@ -291,7 +291,8 @@ namespace {
}
std::string wellhead_location(const context&, std::size_t, std::size_t) const {
auto i { std::to_string(well.getHeadI() + 1) }, j { std::to_string(well.getHeadJ() + 1) } ;
auto i = std::to_string(well.getHeadI() + 1);
auto j = std::to_string(well.getHeadJ() + 1);
right_align(i, 3);
right_align(j, 3);
@ -300,7 +301,10 @@ namespace {
}
std::string reference_depth(const context& ctx, std::size_t, std::size_t) const {
return format_number(ctx.unit_system, Opm::UnitSystem::measure::length, well.getRefDepth(), 6);
if (well.hasRefDepth())
return format_number(ctx.unit_system, Opm::UnitSystem::measure::length, well.getRefDepth(), 6);
else
return format_number(ctx.unit_system, Opm::UnitSystem::measure::identity, -1.0e+20, 9);
}
std::string preferred_phase(const context&, std::size_t, std::size_t) const {
@ -327,10 +331,8 @@ namespace {
return well.segmented_density_calculation() ? "SEG" : "AVG";
}
/*
Don't know what the D-FACTOR represents, but all examples just show 0;
we have therefor hardcoded that for now.
*/
// Don't know what the D-FACTOR represents, but all examples just
// show 0; we have therefore hardcoded that for now.
std::string D_factor(const context&, std::size_t, std::size_t) const {
return "0";
}
@ -363,8 +365,8 @@ namespace {
{ 3, { "PVT" , "TAB" , }, &WellWrapper::pvt_tab , },
{ 4, { "WELL" , "DENS" , "CALC" }, &WellWrapper::dens_calc , },
{ 3, { "FIP" , "REG" , }, &WellWrapper::region_number , },
{ 11, { "WELL" , "D-FACTOR 1" , "DAY/SM3" }, &WellWrapper::D_factor , }};
{ 11, { "WELL" , "D-FACTOR 1" , "DAY/SM3" }, &WellWrapper::D_factor , },
};
void report_well_specification_data(std::ostream& os, const std::vector<Opm::Well>& data, const context& ctx) {

View File

@ -177,7 +177,7 @@ Well::Well(const RestartIO::RstWell& rst_well,
init_step(report_step),
headI(rst_well.ij[0]),
headJ(rst_well.ij[1]),
ref_depth(rst_well.datum_depth),
ref_depth((std::abs(rst_well.datum_depth) < 1.0e+20) ? std::optional<double>{ rst_well.datum_depth } : std::nullopt),
drainage_radius(rst_well.drainage_radius),
allow_cross_flow(rst_well.allow_xflow == 1),
automatic_shutin(def_automatic_shutin),
@ -847,8 +847,13 @@ bool Well::getAllowCrossFlow() const {
return this->allow_cross_flow;
}
bool Well::hasRefDepth() const
{
return this->ref_depth.has_value();
}
double Well::getRefDepth() const {
if (!this->ref_depth.has_value())
if (!this->hasRefDepth())
throw std::logic_error(fmt::format("Well: {} - tried to access not initialized well reference depth", this->name()));
return *this->ref_depth;
}
@ -1621,48 +1626,53 @@ Well::GuideRateTarget Well::GuideRateTargetFromString( const std::string& string
bool Well::cmp_structure(const Well& other) const {
if ((segments && !other.segments) || (!segments && other.segments)) {
if ((this->segments && !other.segments) ||
(!this->segments && other.segments))
{
return false;
}
if (segments && (this->getSegments() != other.getSegments())) {
if (this->segments && (this->getSegments() != other.getSegments())) {
return false;
}
return this->name() == other.name() &&
this->groupName() == other.groupName() &&
this->firstTimeStep() == other.firstTimeStep() &&
this->seqIndex() == other.seqIndex() &&
this->getHeadI() == other.getHeadI() &&
this->getHeadJ() == other.getHeadJ() &&
this->getRefDepth() == other.getRefDepth() &&
this->getPreferredPhase() == other.getPreferredPhase() &&
this->unit_system == other.unit_system &&
this->udq_undefined == other.udq_undefined &&
this->getConnections() == other.getConnections() &&
this->getDrainageRadius() == other.getDrainageRadius() &&
this->getAllowCrossFlow() == other.getAllowCrossFlow() &&
this->getAutomaticShutIn() == other.getAutomaticShutIn() &&
this->getEfficiencyFactor() == other.getEfficiencyFactor();
return (this->name() == other.name())
&& (this->groupName() == other.groupName())
&& (this->firstTimeStep() == other.firstTimeStep())
&& (this->seqIndex() == other.seqIndex())
&& (this->getHeadI() == other.getHeadI())
&& (this->getHeadJ() == other.getHeadJ())
&& (this->hasRefDepth() == other.hasRefDepth())
&& (!this->hasRefDepth() || (this->getRefDepth() == other.getRefDepth()))
&& (this->getPreferredPhase() == other.getPreferredPhase())
&& (this->unit_system == other.unit_system)
&& (this->udq_undefined == other.udq_undefined)
&& (this->getConnections() == other.getConnections())
&& (this->getDrainageRadius() == other.getDrainageRadius())
&& (this->getAllowCrossFlow() == other.getAllowCrossFlow())
&& (this->getAutomaticShutIn() == other.getAutomaticShutIn())
&& (this->getEfficiencyFactor() == other.getEfficiencyFactor())
;
}
bool Well::operator==(const Well& data) const {
return this->cmp_structure(data) &&
this->getSolventFraction() == data.getSolventFraction() &&
this->getEconLimits() == data.getEconLimits() &&
this->isProducer() == data.isProducer() &&
this->getFoamProperties() == data.getFoamProperties() &&
this->getStatus() == data.getStatus() &&
this->guide_rate == data.guide_rate &&
this->solvent_fraction == data.solvent_fraction &&
this->hasProduced() == data.hasProduced() &&
this->hasInjected() == data.hasInjected() &&
this->predictionMode() == data.predictionMode() &&
this->getTracerProperties() == data.getTracerProperties() &&
this->getProductionProperties() == data.getProductionProperties() &&
this->m_pavg == data.m_pavg &&
this->getInjectionProperties() == data.getInjectionProperties();
return this->cmp_structure(data)
&& (this->getSolventFraction() == data.getSolventFraction())
&& (this->getEconLimits() == data.getEconLimits())
&& (this->isProducer() == data.isProducer())
&& (this->getFoamProperties() == data.getFoamProperties())
&& (this->getStatus() == data.getStatus())
&& (this->guide_rate == data.guide_rate)
&& (this->solvent_fraction == data.solvent_fraction)
&& (this->hasProduced() == data.hasProduced())
&& (this->hasInjected() == data.hasInjected())
&& (this->predictionMode() == data.predictionMode())
&& (this->getTracerProperties() == data.getTracerProperties())
&& (this->getProductionProperties() == data.getProductionProperties())
&& (this->m_pavg == data.m_pavg)
&& (this->getInjectionProperties() == data.getInjectionProperties())
;
}

View File

@ -1578,3 +1578,76 @@ END
BOOST_CHECK_EQUAL(w5.getRefDepth(), grid.getCellDepth(0,0,0));
BOOST_CHECK_EQUAL(w5.getWPaveRefDepth(), 0);
}
BOOST_AUTO_TEST_CASE(Missing_RefDepth) {
const auto deck = Parser{}.parseString(R"(RUNSPEC
START
17 DEC 2021 /
DIMENS
10 10 4 /
GRID
DXV
10*100.0 /
DYV
10*100.0 /
DZV
4*10.0 /
DEPTHZ
121*2000.0 /
PERMX
400*100.0 /
PERMY
400*100.0 /
PERMZ
400*10.0 /
PORO
400*0.3 /
-- Deactivate Cells (1,1,3) And (1,1,4)
ACTNUM
1 99*1
1 99*1
0 99*1
0 99*1
/
SCHEDULE
WELSPECS
'W1' 'G' 1 1 1* 'OIL' 2* 'STOP' 4* /
/
COMPDAT
'W1' 1 1 4 4 'OPEN' 1* 34.720 0.216 3095.832 2* 'Y' 12.828 /
'W1' 1 1 3 3 'OPEN' 1* 34.720 0.216 3095.832 2* 'Y' 12.828 /
/
TSTEP
1 /
COMPDAT
'W1' 1 1 2 2 'OPEN' 1* 25.620 0.216 2086.842 2* 'Y' 8.486 /
/
TSTEP
1 /
END
)");
const auto es = EclipseState{ deck };
const auto sched = Schedule{ deck, es };
const auto& w0 = sched[0].wells("W1");
BOOST_CHECK_MESSAGE(! w0.hasRefDepth(),
R"(Well "W1" must NOT have a BHP reference depth at report=1)");
BOOST_CHECK_THROW(w0.getRefDepth(), std::logic_error);
const auto& w1 = sched[1].wells("W1");
BOOST_CHECK_MESSAGE(w1.hasRefDepth(),
R"(Well "W1" must have a BHP reference depth at report=2)");
BOOST_CHECK_CLOSE(w1.getRefDepth(), es.getInputGrid().getCellDepth(0, 0, 1), 1.0e-8);
}