OilfieldToolsNet/opm-simulators
4
0
Fork 0
mirror of https://github.com/OPM/opm-simulators.git synced 2024-11-24 10:10:18 -06:00
Code Issues Packages Projects Releases Wiki Activity

Merge pull request #3837 from totto82/wvfpexp

Browse Source 
Add option for explicit vfp lookup for problematic wells
This commit is contained in:
Tor Harald Sandve 2022-09-09 10:00:02 +02:00 committed by GitHub
commit f36ac3eeed
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
21 changed files with 155 additions and 45 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

19
examples/printvfp.cpp
View File

@ -29,9 +29,11 @@
#include <opm/simulators/wells/VFPProperties.hpp>
#include <opm/simulators/wells/VFPInjProperties.hpp>
#include <opm/simulators/wells/VFPProdProperties.hpp>
#include <opm/simulators/wells/WellState.hpp>
#include <opm/common/utility/TimeService.hpp>
#include <opm/simulators/wells/VFPHelpers.hpp>
#include <opm/core/props/phaseUsageFromDeck.hpp>
#include <iostream>
#include <iomanip>
@ -51,16 +53,17 @@ struct Setup
Parser parser;
auto deck = parser.parseFile(file);
ecl_state.reset(new EclipseState(deck) );
{
const TableManager table( deck );
const Runspec runspec(deck);
python = std::make_shared<Python>();
schedule.reset( new Schedule(deck, *ecl_state, python));
summary_state.reset( new SummaryState(TimeService::from_time_t(schedule->getStartTime())));
}
const TableManager table( deck );
const Runspec runspec(deck);
python = std::make_shared<Python>();
schedule.reset( new Schedule(deck, *ecl_state, python));
summary_state.reset( new SummaryState(TimeService::from_time_t(schedule->getStartTime())));
const int step = 0;
const auto& sched_state = schedule->operator[](step);
vfp_properties = std::make_unique<VFPProperties>(sched_state.vfpinj(), sched_state.vfpprod());
WellState well_state(phaseUsage(runspec.phases()));
vfp_properties = std::make_unique<VFPProperties>(sched_state.vfpinj(), sched_state.vfpprod(), well_state);
};
};

2
opm/simulators/wells/BlackoilWellModel.hpp
View File

@ -142,7 +142,7 @@ namespace Opm {
BlackoilWellModel(Simulator& ebosSimulator);
void init();
void initWellContainer() override;
void initWellContainer(const int reportStepIdx) override;
/////////////
// <eWoms auxiliary module stuff>

2
opm/simulators/wells/BlackoilWellModelGeneric.cpp
View File

@ -2358,7 +2358,7 @@ runWellPIScaling(const int timeStepIdx,
this->createWellContainer(timeStepIdx);
this->inferLocalShutWells();
this->initWellContainer();
this->initWellContainer(timeStepIdx);
this->calculateProductivityIndexValues(local_deferredLogger);
this->calculateProductivityIndexValuesShutWells(timeStepIdx, local_deferredLogger);

2
opm/simulators/wells/BlackoilWellModelGeneric.hpp
View File

@ -400,7 +400,7 @@ protected:
// create the well container
virtual void createWellContainer(const int time_step) = 0;
virtual void initWellContainer() = 0;
virtual void initWellContainer(const int reportStepIdx) = 0;
virtual void calculateProductivityIndexValuesShutWells(const int reportStepIdx,
DeferredLogger& deferred_logger) = 0;

16
opm/simulators/wells/BlackoilWellModel_impl.hpp
View File

@ -93,11 +93,15 @@ namespace Opm {
template<typename TypeTag>
void
BlackoilWellModel<TypeTag>::
initWellContainer()
initWellContainer(const int reportStepIdx)
{
const uint64_t effective_events_mask = ScheduleEvents::WELL_STATUS_CHANGE
+ ScheduleEvents::NEW_WELL;
const auto& events = schedule()[reportStepIdx].wellgroup_events();
for (auto& wellPtr : this->well_container_) {
const bool well_opened_this_step = report_step_starts_ && events.hasEvent(wellPtr->name(), effective_events_mask);
wellPtr->init(&this->phase_usage_, this->depth_, this->gravity_,
this->local_num_cells_, this->B_avg_);
this->local_num_cells_, this->B_avg_, well_opened_this_step);
}
}
@ -222,7 +226,7 @@ namespace Opm {
{
const auto& sched_state = this->schedule()[timeStepIdx];
// update VFP properties
vfp_properties_.reset(new VFPProperties( sched_state.vfpinj(), sched_state.vfpprod()) );
vfp_properties_.reset(new VFPProperties( sched_state.vfpinj(), sched_state.vfpprod(), this->prevWellState()));
this->initializeWellProdIndCalculators();
if (sched_state.events().hasEvent(ScheduleEvents::Events::WELL_PRODUCTIVITY_INDEX)) {
this->runWellPIScaling(timeStepIdx, local_deferredLogger);
@ -263,7 +267,7 @@ namespace Opm {
// do the initialization for all the wells
// TODO: to see whether we can postpone of the intialization of the well containers to
// optimize the usage of the following several member variables
this->initWellContainer();
this->initWellContainer(reportStepIdx);
// update the updated cell flag
std::fill(is_cell_perforated_.begin(), is_cell_perforated_.end(), false);
@ -390,7 +394,7 @@ namespace Opm {
WellInterfacePtr well = createWellForWellTest(well_name, timeStepIdx, deferred_logger);
// some preparation before the well can be used
well->init(&phase_usage_, depth_, gravity_, local_num_cells_, B_avg_);
well->init(&phase_usage_, depth_, gravity_, local_num_cells_, B_avg_, true);
double well_efficiency_factor = wellEcl.getEfficiencyFactor();
WellGroupHelpers::accumulateGroupEfficiencyFactor(schedule().getGroup(wellEcl.groupName(), timeStepIdx),
@ -1577,7 +1581,7 @@ namespace Opm {
<StandardWell<TypeTag>>(shutWell, reportStepIdx);
wellPtr->init(&this->phase_usage_, this->depth_, this->gravity_,
this->local_num_cells_, this->B_avg_);
this->local_num_cells_, this->B_avg_, true);
this->calculateProductivityIndexValues(wellPtr.get(), deferred_logger);
}

3
opm/simulators/wells/MultisegmentWell.hpp
View File

@ -91,7 +91,8 @@ namespace Opm
const std::vector<double>& depth_arg,
const double gravity_arg,
const int num_cells,
const std::vector< Scalar >& B_avg) override;
const std::vector< Scalar >& B_avg,
const bool changed_to_open_this_step) override;
virtual void initPrimaryVariablesEvaluation() const override;

5
opm/simulators/wells/MultisegmentWell_impl.hpp
View File

@ -94,9 +94,10 @@ namespace Opm
const std::vector<double>& depth_arg,
const double gravity_arg,
const int num_cells,
const std::vector< Scalar >& B_avg)
const std::vector< Scalar >& B_avg,
const bool changed_to_open_this_step)
{
Base::init(phase_usage_arg, depth_arg, gravity_arg, num_cells, B_avg);
Base::init(phase_usage_arg, depth_arg, gravity_arg, num_cells, B_avg, changed_to_open_this_step);
// TODO: for StandardWell, we need to update the perf depth here using depth_arg.
// for MultisegmentWell, it is much more complicated.

3
opm/simulators/wells/StandardWell.hpp
View File

@ -139,7 +139,8 @@ namespace Opm
const std::vector<double>& depth_arg,
const double gravity_arg,
const int num_cells,
const std::vector< Scalar >& B_avg) override;
const std::vector< Scalar >& B_avg,
const bool changed_to_open_this_step) override;
virtual void initPrimaryVariablesEvaluation() const override;

5
opm/simulators/wells/StandardWell_impl.hpp
View File

@ -61,9 +61,10 @@ namespace Opm
const std::vector<double>& depth_arg,
const double gravity_arg,
const int num_cells,
const std::vector< Scalar >& B_avg)
const std::vector< Scalar >& B_avg,
const bool changed_to_open_this_step)
{
Base::init(phase_usage_arg, depth_arg, gravity_arg, num_cells, B_avg);
Base::init(phase_usage_arg, depth_arg, gravity_arg, num_cells, B_avg, changed_to_open_this_step);
this->StdWellEval::init(this->perf_depth_, depth_arg, num_cells, Base::has_polymermw);
}

9
opm/simulators/wells/VFPHelpers.cpp
View File

@ -345,12 +345,19 @@ VFPEvaluation bhp(const VFPProdTable& table,
const double liquid,
const double vapour,
const double thp,
const double alq)
const double alq,
const double explicit_wfr,
const double explicit_gfr,
const bool use_vfpexplicit)
{
//Find interpolation variables
double flo = detail::getFlo(table, aqua, liquid, vapour);
double wfr = detail::getWFR(table, aqua, liquid, vapour);
double gfr = detail::getGFR(table, aqua, liquid, vapour);
if (use_vfpexplicit) {
wfr = explicit_wfr;
gfr = explicit_gfr;
}
//First, find the values to interpolate between
//Recall that flo is negative in Opm, so switch sign.

5
opm/simulators/wells/VFPHelpers.hpp
View File

@ -136,7 +136,10 @@ VFPEvaluation bhp(const VFPProdTable& table,
const double liquid,
const double vapour,
const double thp,
const double alq);
const double alq,
const double explicit_wfr,
const double explicit_gfr,
const bool use_vfpexplicit);
VFPEvaluation bhp(const VFPInjTable& table,
const double aqua,

18
opm/simulators/wells/VFPProdProperties.cpp
View File

@ -87,10 +87,13 @@ double VFPProdProperties::bhp(int table_id,
const double& liquid,
const double& vapour,
const double& thp_arg,
const double& alq) const {
const double& alq,
const double& explicit_wfr,
const double& explicit_gfr,
const bool use_expvfp) const {
const VFPProdTable& table = detail::getTable(m_tables, table_id);
detail::VFPEvaluation retval = detail::bhp(table, aqua, liquid, vapour, thp_arg, alq);
detail::VFPEvaluation retval = detail::bhp(table, aqua, liquid, vapour, thp_arg, alq, explicit_wfr,explicit_gfr, use_expvfp);
return retval.value;
}
@ -145,7 +148,10 @@ EvalWell VFPProdProperties::bhp(const int table_id,
const EvalWell& liquid,
const EvalWell& vapour,
const double& thp,
const double& alq) const
const double& alq,
const double& explicit_wfr,
const double& explicit_gfr,
const bool use_expvfp) const
{
//Get the table
const VFPProdTable& table = detail::getTable(m_tables, table_id);
@ -155,6 +161,10 @@ EvalWell VFPProdProperties::bhp(const int table_id,
EvalWell flo = detail::getFlo(table, aqua, liquid, vapour);
EvalWell wfr = detail::getWFR(table, aqua, liquid, vapour);
EvalWell gfr = detail::getGFR(table, aqua, liquid, vapour);
if (use_expvfp) {
wfr = explicit_wfr;
gfr = explicit_gfr;
}
//First, find the values to interpolate between
//Value of FLO is negative in OPM for producers, but positive in VFP table
@ -175,7 +185,7 @@ EvalWell VFPProdProperties::bhp(const int table_id,
#define INSTANCE(...) \
template __VA_ARGS__ VFPProdProperties::bhp<__VA_ARGS__>(const int, \
const __VA_ARGS__&, const __VA_ARGS__&, const __VA_ARGS__&, \
const double&, const double&) const;
const double&, const double&, const double&, const double&, const bool) const;
INSTANCE(DenseAd::Evaluation<double, -1, 4u>)
INSTANCE(DenseAd::Evaluation<double, -1, 5u>)

10
opm/simulators/wells/VFPProdProperties.hpp
View File

@ -65,7 +65,10 @@ public:
const EvalWell& liquid,
const EvalWell& vapour,
const double& thp,
const double& alq) const;
const double& alq,
const double& explicit_wfr,
const double& explicit_gfr,
const bool use_expvfp) const;
/**
* Linear interpolation of bhp as a function of the input parameters
@ -84,7 +87,10 @@ public:
const double& liquid,
const double& vapour,
const double& thp,
const double& alq) const;
const double& alq,
const double& explicit_wfr,
const double& explicit_gfr,
const bool use_expvfp) const;
/**
* Linear interpolation of thp as a function of the input parameters

30
opm/simulators/wells/VFPProperties.hpp
View File

@ -22,6 +22,8 @@
#include <opm/simulators/wells/VFPInjProperties.hpp>
#include <opm/simulators/wells/VFPProdProperties.hpp>
#include <opm/simulators/wells/WellState.hpp>
#include <opm/simulators/wells/VFPHelpers.hpp>
#include <map>
@ -47,7 +49,9 @@ public:
*/
VFPProperties(const std::vector<std::reference_wrapper<const VFPInjTable>>& inj_tables,
const std::vector<std::reference_wrapper<const VFPProdTable>>& prod_tables)
const std::vector<std::reference_wrapper<const VFPProdTable>>& prod_tables,
const WellState& well_state)
:well_state_(well_state)
{
for (const auto& vfpinj : inj_tables)
this->m_inj.addTable( vfpinj );
@ -70,9 +74,33 @@ public:
return &m_prod;
}
double getExplicitWFR(const int table_id, const size_t well_index) const {
const auto& rates = well_state_.well(well_index).surface_rates;
assert(rates.size() == 3);
const auto& pu = well_state_.phaseUsage();
const auto& aqua = rates[pu.phase_pos[BlackoilPhases::Aqua]];
const auto& liquid = rates[pu.phase_pos[BlackoilPhases::Liquid]];
const auto& vapour = rates[pu.phase_pos[BlackoilPhases::Vapour]];
const VFPProdTable& table = this->m_prod.getTable(table_id);
return detail::getWFR(table, aqua, liquid, vapour);
}
double getExplicitGFR(const int table_id, const size_t well_index) const {
const auto& rates = well_state_.well(well_index).surface_rates;
assert(rates.size() == 3);
const auto& pu = well_state_.phaseUsage();
const auto& aqua = rates[pu.phase_pos[BlackoilPhases::Aqua]];
const auto& liquid = rates[pu.phase_pos[BlackoilPhases::Liquid]];
const auto& vapour = rates[pu.phase_pos[BlackoilPhases::Vapour]];
const VFPProdTable& table = this->m_prod.getTable(table_id);
return detail::getGFR(table, aqua, liquid, vapour);
}
private:
VFPInjProperties m_inj;
VFPProdProperties m_prod;
const WellState& well_state_;
};

5
opm/simulators/wells/WellGroupHelpers.cpp
View File

@ -856,7 +856,10 @@ namespace WellGroupHelpers
rates[BlackoilPhases::Liquid],
rates[BlackoilPhases::Vapour],
up_press,
alq);
alq,
0.0, //explicit_wfr
0.0, //explicit_gfr
false); //use_expvfp we dont support explicit lookup
#define EXTRA_DEBUG_NETWORK 0
#if EXTRA_DEBUG_NETWORK
std::ostringstream oss;

3
opm/simulators/wells/WellInterface.hpp
View File

@ -147,7 +147,8 @@ public:
const std::vector<double>& depth_arg,
const double gravity_arg,
const int num_cells,
const std::vector< Scalar >& B_avg);
const std::vector< Scalar >& B_avg,
const bool changed_to_open_this_step);
virtual void initPrimaryVariablesEvaluation() const = 0;

5
opm/simulators/wells/WellInterfaceEval.cpp
View File

@ -493,7 +493,10 @@ calculateBhpFromThp(const WellState& well_state,
const auto& controls = well.productionControls(summaryState);
const double vfp_ref_depth = baseif_.vfpProperties()->getProd()->getTable(controls.vfp_table_number).getDatumDepth();
const double dp = wellhelpers::computeHydrostaticCorrection(baseif_.refDepth(), vfp_ref_depth, rho, baseif_.gravity());
return baseif_.vfpProperties()->getProd()->bhp(controls.vfp_table_number, aqua, liquid, vapour, baseif_.getTHPConstraint(summaryState), baseif_.getALQ(well_state)) - dp;
const auto& wfr = baseif_.vfpProperties()->getExplicitWFR(controls.vfp_table_number, baseif_.indexOfWell());
const auto& gfr = baseif_.vfpProperties()->getExplicitGFR(controls.vfp_table_number, baseif_.indexOfWell());
const bool use_vfpexplicit = baseif_.useVfpExplicit();
return baseif_.vfpProperties()->getProd()->bhp(controls.vfp_table_number, aqua, liquid, vapour, baseif_.getTHPConstraint(summaryState), baseif_.getALQ(well_state), wfr, gfr, use_vfpexplicit) - dp;
}
else {
OPM_DEFLOG_THROW(std::logic_error, "Expected INJECTOR or PRODUCER for well " + baseif_.name(), deferred_logger);

12
opm/simulators/wells/WellInterfaceGeneric.cpp
View File

@ -414,6 +414,12 @@ bool WellInterfaceGeneric::isOperableAndSolvable() const
return operability_status_.isOperableAndSolvable();
}
bool WellInterfaceGeneric::useVfpExplicit() const
{
const auto& wvfpexp = well_ecl_.getWVFPEXP();
return ((wvfpexp.explicit_lookup() && !changedToOpenThisStep())|| operability_status_.use_vfpexplicit);
}
double WellInterfaceGeneric::getALQ(const WellState& well_state) const
{
return well_state.getALQ(name());
@ -656,10 +662,14 @@ computeBhpAtThpLimitProdCommon(const std::function<std::vector<double>(const dou
const double vfp_ref_depth = table.getDatumDepth();
const double thp_limit = this->getTHPConstraint(summary_state);
const double dp = wellhelpers::computeHydrostaticCorrection(this->refDepth(), vfp_ref_depth, rho, this->gravity());
auto fbhp = [this, &controls, thp_limit, dp, alq_value](const std::vector<double>& rates) {
assert(rates.size() == 3);
const auto& wfr = this->vfpProperties()->getExplicitWFR(controls.vfp_table_number, this->indexOfWell());
const auto& gfr = this->vfpProperties()->getExplicitGFR(controls.vfp_table_number, this->indexOfWell());
const bool use_vfpexp = this->useVfpExplicit();
return this->vfpProperties()->getProd()
->bhp(controls.vfp_table_number, rates[Water], rates[Oil], rates[Gas], thp_limit, alq_value) - dp;
->bhp(controls.vfp_table_number, rates[Water], rates[Oil], rates[Gas], thp_limit, alq_value, wfr, gfr, use_vfpexp) - dp;
};
// Make the flo() function.

5
opm/simulators/wells/WellInterfaceGeneric.hpp
View File

@ -88,6 +88,7 @@ public:
// whether the well is operable
bool isOperableAndSolvable() const;
bool useVfpExplicit () const;
void initCompletions();
void closeCompletions(const WellTestState& wellTestState);
@ -265,6 +266,8 @@ protected:
bool has_negative_potentials = false;
//thp limit violated but not switched
mutable bool thp_limit_violated_but_not_switched = false;
bool use_vfpexplicit = false;
};
OperabilityStatus operability_status_;
@ -353,7 +356,7 @@ protected:
std::vector< std::string> well_control_log_;
bool changed_to_open_this_step_ = false;
bool changed_to_open_this_step_ = true;
};
}

12
opm/simulators/wells/WellInterface_impl.hpp
View File

@ -73,11 +73,13 @@ namespace Opm
const std::vector<double>& /* depth_arg */,
const double gravity_arg,
const int /* num_cells */,
const std::vector< Scalar >& B_avg)
const std::vector< Scalar >& B_avg,
const bool changed_to_open_this_step)
{
this->phase_usage_ = phase_usage_arg;
this->gravity_ = gravity_arg;
B_avg_ = B_avg;
this->changed_to_open_this_step_ = changed_to_open_this_step;
}
@ -613,6 +615,12 @@ namespace Opm
}
updateWellOperability(ebos_simulator, well_state, deferred_logger);
if (!this->operability_status_.isOperableAndSolvable()) {
this->operability_status_.use_vfpexplicit = true;
deferred_logger.debug("EXPLICIT_LOOKUP_VFP",
"well not operable, trying with explicit vfp lookup: " + this->name());
updateWellOperability(ebos_simulator, well_state, deferred_logger);
}
}
template<typename TypeTag>
@ -681,7 +689,7 @@ namespace Opm
checkOperabilityUnderBHPLimit(well_state, ebos_simulator, deferred_logger);
}
// we do some extra checking for wells under THP control.
if (thp_controled) {
if (check_thp) {
checkOperabilityUnderTHPLimit(ebos_simulator, well_state, deferred_logger);
}
}

29
tests/test_vfpproperties.cpp
View File

@ -270,7 +270,7 @@ BOOST_AUTO_TEST_CASE(InterpolateZero)
const double v = m / static_cast<double>(n-1);
//Note order of arguments!
sum += properties->bhp(1, v, x, y, z, u);
sum += properties->bhp(1, v, x, y, z, u, 0.0, 0.0, false);
}
}
}
@ -305,7 +305,7 @@ BOOST_AUTO_TEST_CASE(InterpolateOne)
const double v = m / static_cast<double>(n-1);
//Note order of arguments!
const double value = properties->bhp(1, v, x, y, z, u);
const double value = properties->bhp(1, v, x, y, z, u, 0, 0, false);
sum += value;
}
@ -527,7 +527,7 @@ BOOST_AUTO_TEST_CASE(THPToBHPAndBackPlane)
double thp = 0.5;
double alq = 32.9;
double bhp_val = properties->bhp(1, aqua, liquid, vapour, thp, alq);
double bhp_val = properties->bhp(1, aqua, liquid, vapour, thp, alq, 0, 0 , false);
double thp_val = properties->thp(1, aqua, liquid, vapour, bhp_val, alq);
BOOST_CHECK_CLOSE(thp_val, thp, max_d_tol);
@ -546,13 +546,30 @@ BOOST_AUTO_TEST_CASE(THPToBHPAndBackNonTrivial)
double thp = 0.5;
double alq = 32.9;
double bhp_val = properties->bhp(1, aqua, liquid, vapour, thp, alq);
double bhp_val = properties->bhp(1, aqua, liquid, vapour, thp, alq, 0, 0, false);
double thp_val = properties->thp(1, aqua, liquid, vapour, bhp_val, alq);
BOOST_CHECK_CLOSE(thp_val, thp, max_d_tol);
}
BOOST_AUTO_TEST_CASE(ExplicitBhpLookup)
{
fillDataRandom();
initProperties();
double aqua = -0.5;
double liquid = -0.9;
double vapour = -0.1;
double thp = 0.5;
double alq = 32.9;
double wfr = Opm::detail::getWFR(*table, aqua, liquid, vapour);
double gfr = Opm::detail::getGFR(*table, aqua, liquid, vapour);
double bhp_val = properties->bhp(1, aqua, liquid, vapour, thp, alq, wfr, gfr, false);
double bhp_val_explicit = properties->bhp(1, 2*aqua, liquid, 3*vapour, thp, alq, wfr, gfr, true);
BOOST_CHECK_CLOSE(bhp_val, bhp_val_explicit, max_d_tol);
}
BOOST_AUTO_TEST_SUITE_END() // Trivial tests
@ -614,7 +631,7 @@ VFPPROD \n\
double thp = t * 456.78;
double alq = a * 42.24;
double bhp_interp = properties.bhp(42, aqua, liquid, vapour, thp, alq);
double bhp_interp = properties.bhp(42, aqua, liquid, vapour, thp, alq, 0, 0, false);
double bhp_ref = thp;
double thp_interp = properties.thp(42, aqua, liquid, vapour, bhp_ref, alq);
double thp_ref = thp;
@ -711,7 +728,7 @@ BOOST_AUTO_TEST_CASE(ParseInterpolateRealisticVFPPROD)
}
else {
//Value given as pascal, convert to barsa for comparison with reference
double value_i = properties.bhp(32, aqua, liquid, vapour, t_i, a_i) * 10.0e-6;
double value_i = properties.bhp(32, aqua, liquid, vapour, t_i, a_i, 0, 0, false) * 10.0e-6;
double abs_diff = std::abs(value_i - reference[i]);
sad += abs_diff;