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

Merge pull request #5308 from akva2/wellinterface_template_scalar

Browse Source 
WellInterface: template Scalar type
This commit is contained in:
Atgeirr Flø Rasmussen 2024-05-14 14:11:45 +02:00 committed by GitHub
commit ae61361061
No known key found for this signature in database
GPG Key ID: B5690EEEBB952194
34 changed files with 628 additions and 557 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

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

@ -1334,14 +1334,14 @@ calculateEfficiencyFactors(const int reportStepIdx)
}
template<class Scalar>
WellInterfaceGeneric*
WellInterfaceGeneric<Scalar>*
BlackoilWellModelGeneric<Scalar>::
getGenWell(const std::string& well_name)
{
// finding the iterator of the well in wells_ecl
auto well = std::find_if(well_container_generic_.begin(),
well_container_generic_.end(),
[&well_name](const WellInterfaceGeneric* elem)->bool {
[&well_name](const WellInterfaceGeneric<Scalar>* elem)->bool {
return elem->name() == well_name;
});

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

@ -64,7 +64,7 @@ namespace Opm {
struct SimulatorUpdate;
class SummaryConfig;
class VFPProperties;
class WellInterfaceGeneric;
template<class Scalar> class WellInterfaceGeneric;
template<class Scalar> class WellState;
} // namespace Opm
@ -84,7 +84,7 @@ class BlackoilWellModelGeneric
public:
// --------- Types ---------
using GLiftOptWells = std::map<std::string, std::unique_ptr<GasLiftSingleWellGeneric<Scalar>>>;
using GLiftProdWells = std::map<std::string, const WellInterfaceGeneric*>;
using GLiftProdWells = std::map<std::string, const WellInterfaceGeneric<Scalar>*>;
using GLiftWellStateMap = std::map<std::string, std::unique_ptr<GasLiftWellState<Scalar>>>;
BlackoilWellModelGeneric(Schedule& schedule,
@ -115,7 +115,7 @@ public:
const Schedule& schedule() const { return schedule_; }
const PhaseUsage& phaseUsage() const { return phase_usage_; }
const GroupState<Scalar>& groupState() const { return this->active_wgstate_.group_state; }
std::vector<const WellInterfaceGeneric*> genericWells() const
std::vector<const WellInterfaceGeneric<Scalar>*> genericWells() const
{ return {well_container_generic_.begin(), well_container_generic_.end()}; }
/*
@ -542,7 +542,7 @@ protected:
std::function<bool(const Well&)> not_on_process_{};
// a vector of all the wells.
std::vector<WellInterfaceGeneric*> well_container_generic_{};
std::vector<WellInterfaceGeneric<Scalar>*> well_container_generic_{};
std::vector<int> local_shut_wells_{};
@ -589,7 +589,7 @@ protected:
std::map<std::string, std::pair<std::string, std::string>> closed_offending_wells_;
private:
WellInterfaceGeneric* getGenWell(const std::string& well_name);
WellInterfaceGeneric<Scalar>* getGenWell(const std::string& well_name);
};

4
opm/simulators/wells/BlackoilWellModelGuideRates.cpp
View File

@ -582,7 +582,7 @@ guideRateUpdateIsNeeded(const int reportStepIdx) const
const auto& genWells = wellModel_.genericWells();
auto need_update =
std::any_of(genWells.begin(), genWells.end(),
[](const WellInterfaceGeneric* well)
[](const WellInterfaceGeneric<Scalar>* well)
{
return well->changedToOpenThisStep();
});
@ -594,7 +594,7 @@ guideRateUpdateIsNeeded(const int reportStepIdx) const
+ ScheduleEvents::NEW_WELL;
need_update = std::any_of(genWells.begin(), genWells.end(),
[&events](const WellInterfaceGeneric* well)
[&events](const WellInterfaceGeneric<Scalar>* well)
{
return events.hasEvent(well->name(), effective_events_mask);
});

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

@ -51,7 +51,7 @@ public:
const Parallel::Communication& comm,
bool glift_debug);
const WellInterfaceGeneric& getWell() const override { return well_; }
const WellInterfaceGeneric<Scalar>& getWell() const override { return well_; }
private:
std::optional<Scalar>

4
opm/simulators/wells/GasLiftSingleWellGeneric.hpp
View File

@ -40,7 +40,7 @@ class GasLiftWell;
template<class Scalar> class GasLiftWellState;
class Schedule;
class SummaryState;
class WellInterfaceGeneric;
template<class Scalar> class WellInterfaceGeneric;
template<class Scalar> class WellState;
template<class Scalar> class GroupState;
@ -107,7 +107,7 @@ public:
std::unique_ptr<GasLiftWellState<Scalar>> runOptimize(const int iteration_idx);
virtual const WellInterfaceGeneric& getWell() const = 0;
virtual const WellInterfaceGeneric<Scalar>& getWell() const = 0;
protected:
GasLiftSingleWellGeneric(DeferredLogger& deferred_logger,

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

@ -238,7 +238,7 @@ checkThpControl_() const
const Well::ProducerCMode& control_mode =
this->well_state_.well(well_index).production_cmode;
bool thp_control = control_mode == Well::ProducerCMode::THP;
const WellInterfaceGeneric &well = getWell();
const auto& well = getWell();
thp_control = thp_control || well.thpLimitViolatedButNotSwitched();
if (this->debug) {
if (!thp_control) {

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

@ -1118,7 +1118,7 @@ computeDelta(const std::string& well_name)
const GradInfo& gi = this->parent.dec_grads_.at(well_name);
GasLiftWellState<Scalar>& state = *(this->parent.well_state_map_.at(well_name).get());
GasLiftSingleWell& gs_well = *(this->parent.stage1_wells_.at(well_name).get());
const WellInterfaceGeneric& well = gs_well.getWell();
const WellInterfaceGeneric<Scalar>& well = gs_well.getWell();
// only get deltas for wells owned by this rank
if (this->parent.well_state_.wellIsOwned(well.indexOfWell(), well_name)) {
const auto& well_ecl = well.wellEcl();

4
opm/simulators/wells/GasLiftStage2.hpp
View File

@ -38,7 +38,7 @@ template<class Scalar> class GasLiftWellState;
class Group;
template<class Scalar> class GroupState;
class Schedule;
class WellInterfaceGeneric;
template<class Scalar> class WellInterfaceGeneric;
template<class Scalar> class WellState;
template<class Scalar>
@ -46,7 +46,7 @@ class GasLiftStage2 : public GasLiftCommon<Scalar>
{
using GasLiftSingleWell = GasLiftSingleWellGeneric<Scalar>;
using GLiftOptWells = std::map<std::string,std::unique_ptr<GasLiftSingleWell>>;
using GLiftProdWells = std::map<std::string,const WellInterfaceGeneric*>;
using GLiftProdWells = std::map<std::string,const WellInterfaceGeneric<Scalar>*>;
using GLiftWellStateMap = std::map<std::string,std::unique_ptr<GasLiftWellState<Scalar>>>;
using GradPair = std::pair<std::string, Scalar>;
using GradPairItr = typename std::vector<GradPair>::iterator;

4
opm/simulators/wells/MultisegmentWellEquations.cpp
View File

@ -311,7 +311,7 @@ extractCPRPressureMatrix(PressureMatrix& jacobian,
const BVector& weights,
const int pressureVarIndex,
const bool /*use_well_weights*/,
const WellInterfaceGeneric& well,
const WellInterfaceGeneric<Scalar>& well,
const int seg_pressure_var_ind,
const WellState<Scalar>& well_state) const
{
@ -395,7 +395,7 @@ template void MultisegmentWellEquations<double,numWellEq,numEq>:: \
const MultisegmentWellEquations<double,numWellEq,numEq>::BVector&, \
const int, \
const bool, \
const WellInterfaceGeneric&, \
const WellInterfaceGeneric<double>&, \
const int, \
const WellState<double>&) const;

4
opm/simulators/wells/MultisegmentWellEquations.hpp
View File

@ -41,7 +41,7 @@ template<class Scalar> class MultisegmentWellGeneric;
#if COMPILE_BDA_BRIDGE
class WellContributions;
#endif
class WellInterfaceGeneric;
template<class Scalar> class WellInterfaceGeneric;
template<class Scalar> class WellState;
template<class Scalar, int numWellEq, int numEq>
@ -118,7 +118,7 @@ public:
const BVector& weights,
const int pressureVarIndex,
const bool /*use_well_weights*/,
const WellInterfaceGeneric& well,
const WellInterfaceGeneric<Scalar>& well,
const int seg_pressure_var_ind,
const WellState<Scalar>& well_state) const;

7
opm/simulators/wells/MultisegmentWellGeneric.cpp
View File

@ -34,18 +34,15 @@
#include <opm/simulators/wells/WellInterfaceGeneric.hpp>
#include <opm/simulators/wells/WellState.hpp>
#include <cassert>
#include <cmath>
#include <stdexcept>
#include <fmt/format.h>
namespace Opm
{
namespace Opm {
template<typename Scalar>
MultisegmentWellGeneric<Scalar>::
MultisegmentWellGeneric(WellInterfaceGeneric& baseif)
MultisegmentWellGeneric(WellInterfaceGeneric<Scalar>& baseif)
: baseif_(baseif)
{
}

6
opm/simulators/wells/MultisegmentWellGeneric.hpp
View File

@ -32,7 +32,7 @@ namespace Opm
class DeferredLogger;
class SummaryState;
class WellInterfaceGeneric;
template<class Scalar> class WellInterfaceGeneric;
enum class WellSegmentCompPressureDrop;
class WellSegments;
template<class Scalar> class WellState;
@ -52,7 +52,7 @@ public:
int numberOfSegments() const;
protected:
MultisegmentWellGeneric(WellInterfaceGeneric& baseif);
MultisegmentWellGeneric(WellInterfaceGeneric<Scalar>& baseif);
// scale the segment rates and pressure based on well rates and bhp
void scaleSegmentRatesWithWellRates(const std::vector<std::vector<int>>& segment_inlets,
@ -75,7 +75,7 @@ protected:
const double density,
const std::vector<double>& seg_dp) const;
const WellInterfaceGeneric& baseif_;
const WellInterfaceGeneric<Scalar>& baseif_;
};
}

3
opm/simulators/wells/MultisegmentWellSegments.cpp
View File

@ -45,7 +45,6 @@
#include <fmt/format.h>
#include <algorithm>
#include <array>
#include <cmath>
#include <cstddef>
@ -61,7 +60,7 @@ namespace Opm
template<class FluidSystem, class Indices>
MultisegmentWellSegments<FluidSystem,Indices>::
MultisegmentWellSegments(const int numSegments,
WellInterfaceGeneric& well)
WellInterfaceGeneric<Scalar>& well)
: perforations_(numSegments)
, perforation_depth_diffs_(well.numPerfs(), 0.0)
, inlets_(well.wellEcl().getSegments().size())

6
opm/simulators/wells/MultisegmentWellSegments.hpp
View File

@ -33,7 +33,7 @@ namespace Opm {
struct PhaseUsage;
template<class Scalar> class SegmentState;
class UnitSystem;
class WellInterfaceGeneric;
template<class Scalar> class WellInterfaceGeneric;
class SummaryState;
} // namespace Opm
@ -49,7 +49,7 @@ class MultisegmentWellSegments
public:
MultisegmentWellSegments(const int numSegments,
WellInterfaceGeneric& well);
WellInterfaceGeneric<Scalar>& well);
void computeFluidProperties(const EvalWell& temperature,
const EvalWell& saltConcentration,
@ -172,7 +172,7 @@ private:
std::vector<std::vector<EvalWell>> phase_fractions_;
std::vector<std::vector<EvalWell>> phase_viscosities_;
WellInterfaceGeneric& well_;
WellInterfaceGeneric<Scalar>& well_;
void copyPhaseDensities(const unsigned phaseIdx,
const std::size_t stride,

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

@ -287,7 +287,7 @@ namespace Opm
DeferredLogger& deferred_logger)
{
const auto [compute_potential, bhp_controlled_well] =
this->WellInterfaceGeneric::computeWellPotentials(well_potentials, well_state);
this->WellInterfaceGeneric<Scalar>::computeWellPotentials(well_potentials, well_state);
if (!compute_potential) {
return;

4
opm/simulators/wells/StandardWellEquations.cpp
View File

@ -291,7 +291,7 @@ extractCPRPressureMatrix(PressureMatrix& jacobian,
const BVector& weights,
const int pressureVarIndex,
const bool use_well_weights,
const WellInterfaceGeneric& well,
const WellInterfaceGeneric<Scalar>& well,
const int bhp_var_index,
const WellState<Scalar>& well_state) const
{
@ -413,7 +413,7 @@ template void StandardWellEquations<double,N>:: \
const typename StandardWellEquations<double,N>::BVector&, \
const int, \
const bool, \
const WellInterfaceGeneric&, \
const WellInterfaceGeneric<double>&, \
const int, \
const WellState<double>&) const;

4
opm/simulators/wells/StandardWellEquations.hpp
View File

@ -39,7 +39,7 @@ template<class Scalar, int numEq> class StandardWellEquationAccess;
#if COMPILE_BDA_BRIDGE
class WellContributions;
#endif
class WellInterfaceGeneric;
template<class Scalar> class WellInterfaceGeneric;
template<class Scalar> class WellState;
template<class Scalar, int numEq>
@ -115,7 +115,7 @@ public:
const BVector& weights,
const int pressureVarIndex,
const bool use_well_weights,
const WellInterfaceGeneric& well,
const WellInterfaceGeneric<Scalar>& well,
const int bhp_var_index,
const WellState<Scalar>& well_state) const;

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

@ -1706,7 +1706,7 @@ namespace Opm
DeferredLogger& deferred_logger) // const
{
const auto [compute_potential, bhp_controlled_well] =
this->WellInterfaceGeneric::computeWellPotentials(well_potentials, well_state);
this->WellInterfaceGeneric<Scalar>::computeWellPotentials(well_potentials, well_state);
if (!compute_potential) {
return;

6
opm/simulators/wells/WellBhpThpCalculator.hpp
View File

@ -35,14 +35,14 @@ namespace Opm
class DeferredLogger;
class SummaryState;
class Well;
class WellInterfaceGeneric;
template<class Scalar> class WellInterfaceGeneric;
template<class Scalar> class WellState;
//! \brief Class for computing BHP limits.
class WellBhpThpCalculator {
public:
//! \brief Constructor sets reference to well.
WellBhpThpCalculator(const WellInterfaceGeneric& well) : well_(well) {}
WellBhpThpCalculator(const WellInterfaceGeneric<double>& well) : well_(well) {}
//! \brief Checks if well has THP constraints.
bool wellHasTHPConstraints(const SummaryState& summaryState) const;
@ -168,7 +168,7 @@ private:
const double dp,
DeferredLogger& deferred_logger) const;
const WellInterfaceGeneric& well_; //!< Reference to well interface
const WellInterfaceGeneric<double>& well_; //!< Reference to well interface
};
}

6
opm/simulators/wells/WellConstraints.hpp
View File

@ -39,7 +39,7 @@ using RegionId = int;
class Rates;
template<class Scalar> class SingleWellState;
class SummaryState;
class WellInterfaceGeneric;
template<class Scalar> class WellInterfaceGeneric;
enum class WellInjectorCMode;
enum class WellProducerCMode;
@ -47,7 +47,7 @@ enum class WellProducerCMode;
class WellConstraints {
public:
//! \brief Constructor sets reference to well.
WellConstraints(const WellInterfaceGeneric& well) : well_(well) {}
WellConstraints(const WellInterfaceGeneric<double>& well) : well_(well) {}
using RateConvFunc = std::function<void(const RegionId, const int,
const std::vector<double>&,
@ -78,7 +78,7 @@ private:
DeferredLogger& deferred_logger,
const std::optional<Well::ProductionControls>& prod_controls = std::nullopt) const;
const WellInterfaceGeneric& well_; //!< Reference to well interface
const WellInterfaceGeneric<double>& well_; //!< Reference to well interface
};
}

6
opm/simulators/wells/WellConvergence.hpp
View File

@ -30,13 +30,13 @@ namespace Opm
class ConvergenceReport;
class DeferredLogger;
class WellInterfaceGeneric;
template<class Scalar> class WellInterfaceGeneric;
template<class Scalar> class WellState;
class WellConvergence
{
public:
WellConvergence(const WellInterfaceGeneric& well)
WellConvergence(const WellInterfaceGeneric<double>& well)
: well_(well)
{}
@ -62,7 +62,7 @@ public:
ConvergenceReport& report) const;
private:
const WellInterfaceGeneric& well_;
const WellInterfaceGeneric<double>& well_;
};
}

4
opm/simulators/wells/WellFilterCake.cpp
View File

@ -36,7 +36,7 @@ namespace Opm {
template<class Scalar>
void WellFilterCake<Scalar>::
updateFiltrationParticleVolume(const WellInterfaceGeneric& well,
updateFiltrationParticleVolume(const WellInterfaceGeneric<Scalar>& well,
const double dt,
const Scalar conc,
const std::size_t water_index,
@ -78,7 +78,7 @@ updateFiltrationParticleVolume(const WellInterfaceGeneric& well,
template<class Scalar>
void WellFilterCake<Scalar>::
updateInjFCMult(const WellInterfaceGeneric& well,
updateInjFCMult(const WellInterfaceGeneric<Scalar>& well,
WellState<Scalar>& well_state,
DeferredLogger& deferred_logger)
{

6
opm/simulators/wells/WellFilterCake.hpp
View File

@ -26,7 +26,7 @@
namespace Opm {
class DeferredLogger;
class WellInterfaceGeneric;
template<class Scalar> class WellInterfaceGeneric;
template<class Scalar> class WellState;
//! \brief Class for well calculations related to filter cakes.
@ -35,14 +35,14 @@ class WellFilterCake {
public:
//! \brief Update the water injection volume.
//! \details Used for calculation related to cake filtration due to injection activity.
void updateFiltrationParticleVolume(const WellInterfaceGeneric& well,
void updateFiltrationParticleVolume(const WellInterfaceGeneric<Scalar>& well,
const double dt,
const Scalar conc,
const std::size_t water_index,
WellState<Scalar>& well_state);
//! \brief Update the multiplier for well transmissbility due to cake filtration.
void updateInjFCMult(const WellInterfaceGeneric& well,
void updateInjFCMult(const WellInterfaceGeneric<Scalar>& well,
WellState<Scalar>& well_state,
DeferredLogger& deferred_logger);

11
opm/simulators/wells/WellGroupConstraints.hpp
View File

@ -40,16 +40,19 @@ enum class InjectorType;
using RegionId = int;
class Schedule;
class SummaryState;
class WellInterfaceGeneric;
template<class Scalar> class WellInterfaceGeneric;
template<class Scalar> class WellState;
//! \brief Class for computing well group constraints.
class WellGroupConstraints {
public:
//! \brief Constructor sets reference to well.
WellGroupConstraints(const WellInterfaceGeneric& well) : well_(well) {}
WellGroupConstraints(const WellInterfaceGeneric<double>& well) : well_(well) {}
using RateConvFunc = std::function<void(const RegionId, const int, const std::optional<std::string>&, std::vector<double>&)>;
using RateConvFunc = std::function<void(const RegionId,
const int,
const std::optional<std::string>&,
std::vector<double>&)>;
bool checkGroupConstraints(WellState<double>& well_state,
const GroupState<double>& group_state,
@ -79,7 +82,7 @@ private:
const RateConvFunc& rateConverter,
DeferredLogger& deferred_logger) const;
const WellInterfaceGeneric& well_; //!< Reference to well interface
const WellInterfaceGeneric<double>& well_; //!< Reference to well interface
};
}

6
opm/simulators/wells/WellGroupControls.hpp
View File

@ -39,14 +39,14 @@ enum class InjectorType;
using RegionId = int;
class Schedule;
class SummaryState;
class WellInterfaceGeneric;
template<class Scalar> class WellInterfaceGeneric;
template<class Scalar> class WellState;
//! \brief Class for computing well group controls.
class WellGroupControls {
public:
//! \brief Constructor sets reference to well.
WellGroupControls(const WellInterfaceGeneric& well) : well_(well) {}
WellGroupControls(const WellInterfaceGeneric<double>& well) : well_(well) {}
using RateConvFunc = std::function<void(const RegionId, const int, const std::optional<std::string>&, std::vector<double>&)>;
@ -98,7 +98,7 @@ public:
DeferredLogger& deferred_logger) const;
private:
const WellInterfaceGeneric& well_; //!< Reference to well interface
const WellInterfaceGeneric<double>& well_; //!< Reference to well interface
};
}

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

@ -98,7 +98,7 @@ public:
using MatrixBlockType = Dune::FieldMatrix<Scalar, Indices::numEq, Indices::numEq>;
using Eval = typename Base::Eval;
using BVector = Dune::BlockVector<VectorBlockType>;
using PressureMatrix = Dune::BCRSMatrix<Opm::MatrixBlock<double, 1, 1>>;
using PressureMatrix = Dune::BCRSMatrix<Opm::MatrixBlock<Scalar, 1, 1>>;
using RateConverterType =
typename WellInterfaceFluidSystem<FluidSystem>::RateConverterType;
@ -148,17 +148,17 @@ public:
virtual ~WellInterface() = default;
virtual void init(const PhaseUsage* phase_usage_arg,
const std::vector<double>& depth_arg,
const double gravity_arg,
const std::vector<Scalar>& depth_arg,
const Scalar 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() = 0;
virtual ConvergenceReport getWellConvergence(const SummaryState& summary_state,
const WellState<Scalar>& well_state,
const std::vector<double>& B_avg,
const std::vector<Scalar>& B_avg,
DeferredLogger& deferred_logger,
const bool relax_tolerance) const = 0;
@ -186,19 +186,16 @@ public:
DeferredLogger& deferred_logger);
virtual void computeWellRatesWithBhp(
const Simulator& simulator,
const double& bhp,
std::vector<double>& well_flux,
DeferredLogger& deferred_logger
) const = 0;
virtual void computeWellRatesWithBhp(const Simulator& ebosSimulator,
const Scalar& bhp,
std::vector<Scalar>& well_flux,
DeferredLogger& deferred_logger) const = 0;
virtual std::optional<double> computeBhpAtThpLimitProdWithAlq(
const Simulator& simulator,
const SummaryState& summary_state,
const double alq_value,
DeferredLogger& deferred_logger
) const = 0;
virtual std::optional<Scalar>
computeBhpAtThpLimitProdWithAlq(const Simulator& ebos_simulator,
const SummaryState& summary_state,
const Scalar alq_value,
DeferredLogger& deferred_logger) const = 0;
/// using the solution x to recover the solution xw for wells and applying
/// xw to update Well State
@ -216,7 +213,7 @@ public:
// TODO: before we decide to put more information under mutable, this function is not const
virtual void computeWellPotentials(const Simulator& simulator,
const WellState<Scalar>& well_state,
std::vector<double>& well_potentials,
std::vector<Scalar>& well_potentials,
DeferredLogger& deferred_logger) = 0;
virtual void updateWellStateWithTarget(const Simulator& simulator,
@ -226,7 +223,7 @@ public:
virtual void computeWellRatesWithBhpIterations(const Simulator& simulator,
const Scalar& bhp,
std::vector<double>& well_flux,
std::vector<Scalar>& well_flux,
DeferredLogger& deferred_logger) const = 0;
bool updateWellStateWithTHPTargetProd(const Simulator& simulator,
@ -245,7 +242,7 @@ public:
const GroupState<Scalar>& group_state,
const Well::InjectionControls& inj_controls,
const Well::ProductionControls& prod_controls,
const double WQTotal,
const Scalar WQTotal,
DeferredLogger& deferred_logger,
const bool fixed_control = false,
const bool fixed_status = false);
@ -263,7 +260,7 @@ public:
WellState<Scalar>& well_state,
DeferredLogger& deferred_logger) const = 0;
virtual double connectionDensity(const int globalConnIdx,
virtual Scalar connectionDensity(const int globalConnIdx,
const int openConnIdx) const = 0;
/// \brief Wether the Jacobian will also have well contributions in it.
@ -298,12 +295,11 @@ public:
const WellState<Scalar>& well_state,
DeferredLogger& deferred_logger);
bool gliftBeginTimeStepWellTestIterateWellEquations(
const Simulator& simulator,
const double dt,
WellState<Scalar>& well_state,
const GroupState<Scalar>& group_state,
DeferredLogger& deferred_logger);
bool gliftBeginTimeStepWellTestIterateWellEquations(const Simulator& ebos_simulator,
const double dt,
WellState<Scalar>& well_state,
const GroupState<Scalar>& group_state,
DeferredLogger& deferred_logger);
void gliftBeginTimeStepWellTestUpdateALQ(const Simulator& simulator,
WellState<Scalar>& well_state,
@ -325,8 +321,9 @@ public:
/// Compute well rates based on current reservoir conditions and well variables.
/// Used in updateWellStateRates().
virtual std::vector<double> computeCurrentWellRates(const Simulator& simulator,
DeferredLogger& deferred_logger) const = 0;
virtual std::vector<Scalar>
computeCurrentWellRates(const Simulator& simulator,
DeferredLogger& deferred_logger) const = 0;
/// Modify the well_state's rates if there is only one nonzero rate.
/// If so, that rate is kept as is, but the others are set proportionally
@ -345,55 +342,50 @@ public:
return connectionRates_;
}
virtual std::vector<double> getPrimaryVars() const
virtual std::vector<Scalar> getPrimaryVars() const
{
return {};
}
virtual int setPrimaryVars(std::vector<double>::const_iterator)
virtual int setPrimaryVars(typename std::vector<Scalar>::const_iterator)
{
return 0;
}
std::vector<double> wellIndex(const int perf,
std::vector<Scalar> wellIndex(const int perf,
const IntensiveQuantities& intQuants,
const double trans_mult,
const SingleWellState<double>& ws) const;
const Scalar trans_mult,
const SingleWellState<Scalar>& ws) const;
void updateConnectionDFactor(const Simulator& simulator,
SingleWellState<double>& ws) const;
SingleWellState<Scalar>& ws) const;
void updateConnectionTransmissibilityFactor(const Simulator& simulator,
SingleWellState<double>& ws) const;
SingleWellState<Scalar>& ws) const;
protected:
// simulation parameters
const ModelParameters& param_;
std::vector<RateVector> connectionRates_;
std::vector< Scalar > B_avg_;
std::vector<Scalar> B_avg_;
bool changed_to_stopped_this_step_ = false;
bool thp_update_iterations = false;
double wpolymer() const;
Scalar wpolymer() const;
Scalar wfoam() const;
Scalar wsalt() const;
Scalar wmicrobes() const;
Scalar woxygen() const;
Scalar wurea() const;
double wfoam() const;
double wsalt() const;
double wmicrobes() const;
double woxygen() const;
double wurea() const;
virtual double getRefDensity() const = 0;
virtual Scalar getRefDensity() const = 0;
// Component fractions for each phase for the well
const std::vector<double>& compFrac() const;
const std::vector<Scalar>& compFrac() const;
std::vector<double> initialWellRateFractions(const Simulator& simulator,
const WellState<Scalar>& well_state) const;
std::vector<Scalar>
initialWellRateFractions(const Simulator& ebosSimulator,
const WellState<Scalar>& well_state) const;
// check whether the well is operable under BHP limit with current reservoir condition
virtual void checkOperabilityUnderBHPLimit(const WellState<Scalar>& well_state,
@ -453,15 +445,16 @@ protected:
const GroupState<Scalar>& group_state,
DeferredLogger& deferred_logger);
std::optional<double> estimateOperableBhp(const Simulator& simulator,
const double dt,
WellState<Scalar>& well_state,
const SummaryState& summary_state,
DeferredLogger& deferred_logger);
std::optional<Scalar>
estimateOperableBhp(const Simulator& ebos_simulator,
const double dt,
WellState<Scalar>& well_state,
const SummaryState& summary_state,
DeferredLogger& deferred_logger);
bool solveWellWithBhp(const Simulator& simulator,
const double dt,
const double bhp,
const Scalar bhp,
WellState<Scalar>& well_state,
DeferredLogger& deferred_logger);
@ -486,21 +479,20 @@ protected:
[[maybe_unused]] DeferredLogger& deferred_logger) const;
void computeConnLevelProdInd(const FluidState& fs,
const std::function<double(const double)>& connPICalc,
const std::function<Scalar(const Scalar)>& connPICalc,
const std::vector<Scalar>& mobility,
double* connPI) const;
Scalar* connPI) const;
void computeConnLevelInjInd(const FluidState& fs,
const Phase preferred_phase,
const std::function<double(const double)>& connIICalc,
const std::function<Scalar(const Scalar)>& connIICalc,
const std::vector<Scalar>& mobility,
double* connII,
Scalar* connII,
DeferredLogger& deferred_logger) const;
double computeConnectionDFactor(const int perf,
Scalar computeConnectionDFactor(const int perf,
const IntensiveQuantities& intQuants,
const SingleWellState<double>& ws) const;
const SingleWellState<Scalar>& ws) const;
};
} // namespace Opm

111
opm/simulators/wells/WellInterfaceFluidSystem.cpp
View File

@ -54,9 +54,9 @@ WellInterfaceFluidSystem(const Well& well,
const int num_phases,
const int index_of_well,
const std::vector<PerforationData>& perf_data)
: WellInterfaceGeneric(well, parallel_well_info, time_step,
pvtRegionIdx, num_components, num_phases,
index_of_well, perf_data)
: WellInterfaceGeneric<Scalar>(well, parallel_well_info, time_step,
pvtRegionIdx, num_components, num_phases,
index_of_well, perf_data)
, rateConverter_(rate_converter)
{
}
@ -64,9 +64,9 @@ WellInterfaceFluidSystem(const Well& well,
template <typename FluidSystem>
void
WellInterfaceFluidSystem<FluidSystem>::
calculateReservoirRates(SingleWellState<double>& ws) const
calculateReservoirRates(SingleWellState<Scalar>& ws) const
{
const int np = number_of_phases_;
const int np = this->number_of_phases_;
const auto& pu = this->phaseUsage();
// Calculate reservoir rates from average pressure and temperature
if ( !pu.has_energy || this->wellEcl().isProducer()) {
@ -82,11 +82,11 @@ calculateReservoirRates(SingleWellState<double>& ws) const
const auto num_perf_well = perf_data.size();
const auto& surf_perf_rates = perf_data.phase_rates;
for (auto i = 0*num_perf_well; i < num_perf_well; ++i) {
const auto surface_rates_perf = std::vector<double>
const auto surface_rates_perf = std::vector<Scalar>
{ surf_perf_rates.begin() + (i + 0)*np ,
surf_perf_rates.begin() + (i + 1)*np };
std::vector<double> voidage_rates_perf(np, 0.0);
std::vector<Scalar> voidage_rates_perf(np, 0.0);
this->rateConverter_
.calcReservoirVoidageRates(fipreg,
this->pvtRegionIdx_,
@ -101,11 +101,11 @@ calculateReservoirRates(SingleWellState<double>& ws) const
}
// For injectors in a thermal case we convert using the well bhp and temperature
// Assume pure phases in the injector
const auto saltConc = 0.0;
auto rsMax = 0.0;
auto rvMax = 0.0;
auto rswMax = 0.0;
auto rvwMax = 0.0;
const Scalar saltConc = 0.0;
Scalar rsMax = 0.0;
Scalar rvMax = 0.0;
Scalar rswMax = 0.0;
Scalar rvwMax = 0.0;
this->rateConverter_
.calcReservoirVoidageRates(this->pvtRegionIdx_,
ws.bhp,
@ -124,14 +124,14 @@ calculateReservoirRates(SingleWellState<double>& ws) const
const auto num_perf_well = perf_data.size();
const auto& surf_perf_rates = perf_data.phase_rates;
for (auto i = 0*num_perf_well; i < num_perf_well; ++i) {
const auto surface_rates_perf = std::vector<double>
const auto surface_rates_perf = std::vector<Scalar>
{ surf_perf_rates.begin() + (i + 0)*np ,
surf_perf_rates.begin() + (i + 1)*np };
const auto pressure = perf_data.pressure[i];
// Calculate other per-phase dynamic quantities.
const auto temperature = ws.temperature; // Assume same temperature in the well
std::vector<double> voidage_rates_perf(np, 0.0);
std::vector<Scalar> voidage_rates_perf(np, 0.0);
this->rateConverter_
.calcReservoirVoidageRates(this->pvtRegionIdx_,
pressure,
@ -153,7 +153,7 @@ calculateReservoirRates(SingleWellState<double>& ws) const
template <typename FluidSystem>
bool
WellInterfaceFluidSystem<FluidSystem>::
checkIndividualConstraints(SingleWellState<double>& ws,
checkIndividualConstraints(SingleWellState<Scalar>& ws,
const SummaryState& summaryState,
DeferredLogger& deferred_logger,
const std::optional<Well::InjectionControls>& inj_controls,
@ -161,8 +161,8 @@ checkIndividualConstraints(SingleWellState<double>& ws,
{
auto rRates = [this](const int fipreg,
const int pvtRegion,
const std::vector<double>& surface_rates,
std::vector<double>& voidage_rates)
const std::vector<Scalar>& surface_rates,
std::vector<Scalar>& voidage_rates)
{
return rateConverter_.calcReservoirVoidageRates(fipreg, pvtRegion,
surface_rates, voidage_rates);
@ -177,23 +177,25 @@ checkIndividualConstraints(SingleWellState<double>& ws,
template <typename FluidSystem>
bool
WellInterfaceFluidSystem<FluidSystem>::
checkGroupConstraints(WellState<double>& well_state,
const GroupState<double>& group_state,
checkGroupConstraints(WellState<Scalar>& well_state,
const GroupState<Scalar>& group_state,
const Schedule& schedule,
const SummaryState& summaryState,
DeferredLogger& deferred_logger) const
{
if (!this->wellEcl().isAvailableForGroupControl())
return false;
auto rCoeff = [this, &group_state](const RegionId id, const int region, const std::optional<std::string>& prod_gname, std::vector<double>& coeff)
auto rCoeff = [this, &group_state](const RegionId id,
const int region,
const std::optional<std::string>& prod_gname,
std::vector<Scalar>& coeff)
{
if (prod_gname)
this->rateConverter().calcCoeff(id, region, group_state.production_rates(*prod_gname), coeff);
this->rateConverter().calcCoeff(id, region,
group_state.production_rates(*prod_gname), coeff);
else
this->rateConverter().calcInjCoeff(id, region, coeff);
};
return WellGroupConstraints(*this).checkGroupConstraints(well_state, group_state,
@ -204,17 +206,19 @@ checkGroupConstraints(WellState<double>& well_state,
template <typename FluidSystem>
bool
WellInterfaceFluidSystem<FluidSystem>::
checkConstraints(WellState<double>& well_state,
const GroupState<double>& group_state,
checkConstraints(WellState<Scalar>& well_state,
const GroupState<Scalar>& group_state,
const Schedule& schedule,
const SummaryState& summaryState,
DeferredLogger& deferred_logger) const
{
const bool ind_broken = checkIndividualConstraints(well_state.well(this->index_of_well_), summaryState, deferred_logger);
const bool ind_broken = checkIndividualConstraints(well_state.well(this->index_of_well_),
summaryState, deferred_logger);
if (ind_broken) {
return true;
} else {
return checkGroupConstraints(well_state, group_state, schedule, summaryState, deferred_logger);
return checkGroupConstraints(well_state, group_state, schedule,
summaryState, deferred_logger);
}
}
@ -236,57 +240,68 @@ flowPhaseToModelPhaseIdx(const int phaseIdx) const
}
template<typename FluidSystem>
std::optional<double>
std::optional<typename WellInterfaceFluidSystem<FluidSystem>::Scalar>
WellInterfaceFluidSystem<FluidSystem>::
getGroupInjectionTargetRate(const Group& group,
const WellState<double>& well_state,
const GroupState<double>& group_state,
const WellState<Scalar>& well_state,
const GroupState<Scalar>& group_state,
const Schedule& schedule,
const SummaryState& summaryState,
const InjectorType& injectorType,
double efficiencyFactor,
Scalar efficiencyFactor,
DeferredLogger& deferred_logger) const
{
auto rCoeff = [this, &group_state](const RegionId id, const int region, const std::optional<std::string>& prod_gname, std::vector<double>& coeff)
auto rCoeff = [this, &group_state](const RegionId id, const int region,
const std::optional<std::string>& prod_gname,
std::vector<Scalar>& coeff)
{
if (prod_gname)
this->rateConverter().calcCoeff(id, region, group_state.production_rates(*prod_gname), coeff);
this->rateConverter().calcCoeff(id, region,
group_state.production_rates(*prod_gname), coeff);
else
this->rateConverter().calcInjCoeff(id, region, coeff);
};
return WellGroupControls(*this).getGroupInjectionTargetRate(group, well_state,
group_state, schedule,
summaryState, injectorType,
rCoeff, efficiencyFactor,
return WellGroupControls(*this).getGroupInjectionTargetRate(group,
well_state,
group_state,
schedule,
summaryState,
injectorType,
rCoeff,
efficiencyFactor,
deferred_logger);
}
template<typename FluidSystem>
double
typename WellInterfaceFluidSystem<FluidSystem>::Scalar
WellInterfaceFluidSystem<FluidSystem>::
getGroupProductionTargetRate(const Group& group,
const WellState<double>& well_state,
const GroupState<double>& group_state,
const WellState<Scalar>& well_state,
const GroupState<Scalar>& group_state,
const Schedule& schedule,
const SummaryState& summaryState,
double efficiencyFactor,
Scalar efficiencyFactor,
DeferredLogger& deferred_logger) const
{
auto rCoeff = [this, &group_state](const RegionId id, const int region, const std::optional<std::string>& prod_gname, std::vector<double>& coeff)
auto rCoeff = [this, &group_state](const RegionId id, const int region,
const std::optional<std::string>& prod_gname,
std::vector<Scalar>& coeff)
{
if (prod_gname)
this->rateConverter().calcCoeff(id, region, group_state.production_rates(*prod_gname), coeff);
this->rateConverter().calcCoeff(id, region,
group_state.production_rates(*prod_gname), coeff);
else
this->rateConverter().calcInjCoeff(id, region, coeff);
};
return WellGroupControls(*this).getGroupProductionTargetRate(group, well_state,
group_state, schedule,
return WellGroupControls(*this).getGroupProductionTargetRate(group,
well_state,
group_state,
schedule,
summaryState,
rCoeff, efficiencyFactor,
rCoeff,
efficiencyFactor,
deferred_logger);
}

33
opm/simulators/wells/WellInterfaceFluidSystem.hpp
View File

@ -44,7 +44,8 @@ template<class Scalar> class SingleWellState;
template<class Scalar> class WellState;
template<class FluidSystem>
class WellInterfaceFluidSystem : public WellInterfaceGeneric {
class WellInterfaceFluidSystem : public WellInterfaceGeneric<typename FluidSystem::Scalar>
{
protected:
using RateConverterType = RateConverter::
SurfaceToReservoirVoidage<FluidSystem, std::vector<int>>;
@ -52,6 +53,8 @@ protected:
static constexpr int INVALIDCOMPLETION = std::numeric_limits<int>::max();
public:
using Scalar = typename FluidSystem::Scalar;
int flowPhaseToModelPhaseIdx(const int phaseIdx) const;
static constexpr int Water = BlackoilPhases::Aqua;
@ -75,43 +78,43 @@ protected:
const std::vector<PerforationData>& perf_data);
// updating the voidage rates in well_state when requested
void calculateReservoirRates(SingleWellState<double>& ws) const;
void calculateReservoirRates(SingleWellState<Scalar>& ws) const;
bool checkIndividualConstraints(SingleWellState<double>& ws,
bool checkIndividualConstraints(SingleWellState<Scalar>& ws,
const SummaryState& summaryState,
DeferredLogger& deferred_logger,
const std::optional<Well::InjectionControls>& inj_controls = std::nullopt,
const std::optional<Well::ProductionControls>& prod_controls = std::nullopt) const;
bool checkGroupConstraints(WellState<double>& well_state,
const GroupState<double>& group_state,
bool checkGroupConstraints(WellState<Scalar>& well_state,
const GroupState<Scalar>& group_state,
const Schedule& schedule,
const SummaryState& summaryState,
DeferredLogger& deferred_logger) const;
bool checkConstraints(WellState<double>& well_state,
const GroupState<double>& group_state,
bool checkConstraints(WellState<Scalar>& well_state,
const GroupState<Scalar>& group_state,
const Schedule& schedule,
const SummaryState& summaryState,
DeferredLogger& deferred_logger) const;
std::optional<double>
std::optional<Scalar>
getGroupInjectionTargetRate(const Group& group,
const WellState<double>& well_state,
const GroupState<double>& group_state,
const WellState<Scalar>& well_state,
const GroupState<Scalar>& group_state,
const Schedule& schedule,
const SummaryState& summaryState,
const InjectorType& injectorType,
double efficiencyFactor,
Scalar efficiencyFactor,
DeferredLogger& deferred_logger) const;
double
Scalar
getGroupProductionTargetRate(const Group& group,
const WellState<double>& well_state,
const GroupState<double>& group_state,
const WellState<Scalar>& well_state,
const GroupState<Scalar>& group_state,
const Schedule& schedule,
const SummaryState& summaryState,
double efficiencyFactor,
Scalar efficiencyFactor,
DeferredLogger& deferred_logger) const;
// For the conversion between the surface volume rate and reservoir voidage rate

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

@ -49,17 +49,18 @@
#include <cstddef>
#include <stdexcept>
namespace Opm
{
namespace Opm {
WellInterfaceGeneric::WellInterfaceGeneric(const Well& well,
const ParallelWellInfo& pw_info,
const int time_step,
const int pvtRegionIdx,
const int num_components,
const int num_phases,
const int index_of_well,
const std::vector<PerforationData>& perf_data)
template<class Scalar>
WellInterfaceGeneric<Scalar>::
WellInterfaceGeneric(const Well& well,
const ParallelWellInfo& pw_info,
const int time_step,
const int pvtRegionIdx,
const int num_components,
const int num_phases,
const int index_of_well,
const std::vector<PerforationData>& perf_data)
: well_ecl_(well)
, parallel_well_info_(pw_info)
, current_step_(time_step)
@ -121,7 +122,9 @@ WellInterfaceGeneric::WellInterfaceGeneric(const Well& well,
// value in VFPPROD record 5 (GFR values) and supplying a dummy input value
// for the gas rate to the methods in VFPProdProperties.cpp, we can extend
// the VFP calculations to the two-phase oil-water case.
void WellInterfaceGeneric::adaptRatesForVFP(std::vector<double>& rates) const
template<class Scalar>
void WellInterfaceGeneric<Scalar>::
adaptRatesForVFP(std::vector<Scalar>& rates) const
{
const auto& pu = this->phaseUsage();
if (pu.num_phases == 2) {
@ -139,73 +142,91 @@ void WellInterfaceGeneric::adaptRatesForVFP(std::vector<double>& rates) const
}
}
const std::vector<PerforationData>& WellInterfaceGeneric::perforationData() const
template<class Scalar>
const std::vector<PerforationData>&
WellInterfaceGeneric<Scalar>::perforationData() const
{
return *perf_data_;
}
const std::string& WellInterfaceGeneric::name() const
template<class Scalar>
const std::string&
WellInterfaceGeneric<Scalar>::name() const
{
return well_ecl_.name();
}
bool WellInterfaceGeneric::isInjector() const
template<class Scalar>
bool WellInterfaceGeneric<Scalar>::isInjector() const
{
return well_ecl_.isInjector();
}
bool WellInterfaceGeneric::isProducer() const
template<class Scalar>
bool WellInterfaceGeneric<Scalar>::isProducer() const
{
return well_ecl_.isProducer();
}
int WellInterfaceGeneric::indexOfWell() const
template<class Scalar>
int WellInterfaceGeneric<Scalar>::indexOfWell() const
{
return index_of_well_;
}
bool WellInterfaceGeneric::getAllowCrossFlow() const
template<class Scalar>
bool WellInterfaceGeneric<Scalar>::getAllowCrossFlow() const
{
return well_ecl_.getAllowCrossFlow();
}
const Well& WellInterfaceGeneric::wellEcl() const
template<class Scalar>
const Well& WellInterfaceGeneric<Scalar>::wellEcl() const
{
return well_ecl_;
}
Well& WellInterfaceGeneric::wellEcl()
template<class Scalar>
Well& WellInterfaceGeneric<Scalar>::wellEcl()
{
return well_ecl_;
}
const PhaseUsage& WellInterfaceGeneric::phaseUsage() const
template<class Scalar>
const PhaseUsage& WellInterfaceGeneric<Scalar>::phaseUsage() const
{
assert(phase_usage_ != nullptr);
return *phase_usage_;
}
double WellInterfaceGeneric::wsolvent() const
template<class Scalar>
Scalar WellInterfaceGeneric<Scalar>::wsolvent() const
{
return wsolvent_;
}
double WellInterfaceGeneric::rsRvInj() const
template<class Scalar>
Scalar WellInterfaceGeneric<Scalar>::rsRvInj() const
{
return well_ecl_.getInjectionProperties().rsRvInj;
}
void WellInterfaceGeneric::initInjMult(const std::vector<double>& max_inj_mult)
template<class Scalar>
void WellInterfaceGeneric<Scalar>::
initInjMult(const std::vector<Scalar>& max_inj_mult)
{
// prev_inj_multiplier_ will stay unchanged during the time step
// while inj_multiplier_ might be updated during the time step
this->prev_inj_multiplier_ = max_inj_mult;
// initializing the inj_multipler_ to be 1.0
this->inj_multiplier_ = std::vector<double>(max_inj_mult.size(), 1.);
this->inj_multiplier_ = std::vector<Scalar>(max_inj_mult.size(), 1.);
}
void WellInterfaceGeneric::updateInjMult(std::vector<double>& inj_multipliers, DeferredLogger& deferred_logger) const
template<class Scalar>
void WellInterfaceGeneric<Scalar>::
updateInjMult(std::vector<Scalar>& inj_multipliers,
DeferredLogger& deferred_logger) const
{
if (inj_multipliers.size() != this->inj_multiplier_.size()) {
OPM_DEFLOG_THROW(std::runtime_error,
@ -217,15 +238,15 @@ void WellInterfaceGeneric::updateInjMult(std::vector<double>& inj_multipliers, D
inj_multipliers = this->inj_multiplier_;
}
double WellInterfaceGeneric::getInjMult(const int perf,
const double bhp,
const double perf_pres) const
template<class Scalar>
Scalar WellInterfaceGeneric<Scalar>::
getInjMult(const int perf,
const Scalar bhp,
const Scalar perf_pres) const
{
assert(!this->isProducer());
double multiplier = 1.;
Scalar multiplier = 1.;
const auto perf_ecl_index = this->perforationData()[perf].ecl_index;
const bool is_wrev = this->well_ecl_.getInjMultMode() == Well::InjMultMode::WREV;
@ -236,7 +257,7 @@ double WellInterfaceGeneric::getInjMult(const int perf,
if (active_injmult) {
const auto& injmult= is_wrev ? this->well_ecl_.getWellInjMult()
: this->well_ecl_.getConnections()[perf_ecl_index].injmult();
const double pres = is_wrev ? bhp : perf_pres;
const Scalar pres = is_wrev ? bhp : perf_pres;
const auto frac_press = injmult.fracture_pressure;
const auto gradient = injmult.multiplier_gradient;
@ -255,10 +276,9 @@ double WellInterfaceGeneric::getInjMult(const int perf,
return multiplier;
}
bool WellInterfaceGeneric::wellHasTHPConstraints(const SummaryState& summaryState) const
template<class Scalar>
bool WellInterfaceGeneric<Scalar>::
wellHasTHPConstraints(const SummaryState& summaryState) const
{
// only wells under prediction mode can have THP constraint
if (!this->wellEcl().predictionMode()) {
@ -272,11 +292,13 @@ bool WellInterfaceGeneric::wellHasTHPConstraints(const SummaryState& summaryStat
return WellBhpThpCalculator(*this).wellHasTHPConstraints(summaryState);
}
void WellInterfaceGeneric::updateWellTestState(const SingleWellState<double>& ws,
const double& simulationTime,
const bool& writeMessageToOPMLog,
WellTestState& wellTestState,
DeferredLogger& deferred_logger) const
template<class Scalar>
void WellInterfaceGeneric<Scalar>::
updateWellTestState(const SingleWellState<Scalar>& ws,
const double& simulationTime,
const bool& writeMessageToOPMLog,
WellTestState& wellTestState,
DeferredLogger& deferred_logger) const
{
// updating well test state based on Economic limits for operable wells
if (this->isOperableAndSolvable()) {
@ -289,7 +311,9 @@ void WellInterfaceGeneric::updateWellTestState(const SingleWellState<double>& ws
// TODO: well can be shut/closed due to other reasons
}
double WellInterfaceGeneric::getTHPConstraint(const SummaryState& summaryState) const
template<class Scalar>
Scalar WellInterfaceGeneric<Scalar>::
getTHPConstraint(const SummaryState& summaryState) const
{
if (dynamic_thp_limit_) {
return *dynamic_thp_limit_;
@ -298,12 +322,14 @@ double WellInterfaceGeneric::getTHPConstraint(const SummaryState& summaryState)
return WellBhpThpCalculator(*this).getTHPConstraint(summaryState);
}
bool WellInterfaceGeneric::underPredictionMode() const
template<class Scalar>
bool WellInterfaceGeneric<Scalar>::underPredictionMode() const
{
return well_ecl_.predictionMode();
}
void WellInterfaceGeneric::initCompletions()
template<class Scalar>
void WellInterfaceGeneric<Scalar>::initCompletions()
{
assert(completions_.empty() );
@ -330,7 +356,9 @@ void WellInterfaceGeneric::initCompletions()
assert(my_next_perf == perf_data_->end());
}
void WellInterfaceGeneric::closeCompletions(const WellTestState& wellTestState)
template<class Scalar>
void WellInterfaceGeneric<Scalar>::
closeCompletions(const WellTestState& wellTestState)
{
const auto& connections = well_ecl_.getConnections();
int perfIdx = 0;
@ -344,46 +372,55 @@ void WellInterfaceGeneric::closeCompletions(const WellTestState& wellTestState)
}
}
void WellInterfaceGeneric::setVFPProperties(const VFPProperties* vfp_properties_arg)
template<class Scalar>
void WellInterfaceGeneric<Scalar>::
setVFPProperties(const VFPProperties* vfp_properties_arg)
{
vfp_properties_ = vfp_properties_arg;
}
void WellInterfaceGeneric::setPrevSurfaceRates(WellState<double>& well_state,
const WellState<double>& prev_well_state) const
{
auto& ws = well_state.well(this->index_of_well_);
auto& ws_prev = prev_well_state.well(this->index_of_well_);
// The logic here is a bit fragile:
// We need non-zero prev_surface_rates for the purpose of providing explicit fractions
// (if needed) for vfp interpolation.
// We assume that current surface rates either are initialized from previous step
// or (if newly opened) from updateWellStateRates. This is fine unless well was
// stopped in previous step in which case it's rates will be zero. In this case,
// we select the previous rates of the previous well state (and hope for the best).
const bool zero_rates = std::all_of(ws.surface_rates.begin(), ws.surface_rates.end(),
[](double rate) {
return rate == 0.; // TODO: should we use a threshhold for comparison?
} );
template<class Scalar>
void WellInterfaceGeneric<Scalar>::
setPrevSurfaceRates(WellState<Scalar>& well_state,
const WellState<Scalar>& prev_well_state) const
{
auto& ws = well_state.well(this->index_of_well_);
auto& ws_prev = prev_well_state.well(this->index_of_well_);
// The logic here is a bit fragile:
// We need non-zero prev_surface_rates for the purpose of providing explicit fractions
// (if needed) for vfp interpolation.
// We assume that current surface rates either are initialized from previous step
// or (if newly opened) from updateWellStateRates. This is fine unless well was
// stopped in previous step in which case it's rates will be zero. In this case,
// we select the previous rates of the previous well state (and hope for the best).
const bool zero_rates = std::all_of(ws.surface_rates.begin(), ws.surface_rates.end(),
[](Scalar rate) {
return rate == 0.; // TODO: should we use a threshhold for comparison?
} );
if (zero_rates) {
ws.prev_surface_rates = ws_prev.prev_surface_rates;
} else {
ws.prev_surface_rates = ws.surface_rates;
}
if (zero_rates) {
ws.prev_surface_rates = ws_prev.prev_surface_rates;
} else {
ws.prev_surface_rates = ws.surface_rates;
}
}
void WellInterfaceGeneric::setGuideRate(const GuideRate* guide_rate_arg)
template<class Scalar>
void WellInterfaceGeneric<Scalar>::
setGuideRate(const GuideRate* guide_rate_arg)
{
guide_rate_ = guide_rate_arg;
}
void WellInterfaceGeneric::setWellEfficiencyFactor(const double efficiency_factor)
template<class Scalar>
void WellInterfaceGeneric<Scalar>::
setWellEfficiencyFactor(const Scalar efficiency_factor)
{
well_efficiency_factor_ = efficiency_factor;
}
void WellInterfaceGeneric::setRepRadiusPerfLength()
template<class Scalar>
void WellInterfaceGeneric<Scalar>::setRepRadiusPerfLength()
{
const int nperf = number_of_perforations_;
@ -411,10 +448,10 @@ void WellInterfaceGeneric::setRepRadiusPerfLength()
assert(my_next_perf->ecl_index == c);
const auto& connection = connections[c];
if (connection.state() == Connection::State::OPEN) {
double radius = connection.rw();
double re = connection.re(); // area equivalent radius of the grid block
double perf_length = connection.connectionLength(); // the length of the well perforation
const double repR = std::sqrt(re * radius);
Scalar radius = connection.rw();
Scalar re = connection.re(); // area equivalent radius of the grid block
Scalar perf_length = connection.connectionLength(); // the length of the well perforation
const Scalar repR = std::sqrt(re * radius);
perf_rep_radius_.push_back(repR);
perf_length_.push_back(perf_length);
bore_diameters_.push_back(2. * radius);
@ -426,30 +463,39 @@ void WellInterfaceGeneric::setRepRadiusPerfLength()
assert(num_active_connections == nperf);
}
void WellInterfaceGeneric::setWsolvent(const double wsolvent)
template<class Scalar>
void WellInterfaceGeneric<Scalar>::
setWsolvent(const Scalar wsolvent)
{
wsolvent_ = wsolvent;
}
void WellInterfaceGeneric::setDynamicThpLimit(const double thp_limit)
template<class Scalar>
void WellInterfaceGeneric<Scalar>::
setDynamicThpLimit(const Scalar thp_limit)
{
dynamic_thp_limit_ = thp_limit;
}
std::optional<double> WellInterfaceGeneric::getDynamicThpLimit() const
template<class Scalar>
std::optional<Scalar>
WellInterfaceGeneric<Scalar>::getDynamicThpLimit() const
{
return dynamic_thp_limit_;
}
void WellInterfaceGeneric::updatePerforatedCell(std::vector<bool>& is_cell_perforated)
template<class Scalar>
void WellInterfaceGeneric<Scalar>::
updatePerforatedCell(std::vector<bool>& is_cell_perforated)
{
for (int perf_idx = 0; perf_idx<number_of_perforations_; ++perf_idx) {
for (int perf_idx = 0; perf_idx < number_of_perforations_; ++perf_idx) {
is_cell_perforated[well_cells_[perf_idx]] = true;
}
}
bool WellInterfaceGeneric::isVFPActive(DeferredLogger& deferred_logger) const
template<class Scalar>
bool WellInterfaceGeneric<Scalar>::
isVFPActive(DeferredLogger& deferred_logger) const
{
// since the well_controls only handles the VFP number when THP constraint/target is there.
// we need to get the table number through the parser, in case THP constraint/target is not there.
@ -489,23 +535,28 @@ bool WellInterfaceGeneric::isVFPActive(DeferredLogger& deferred_logger) const
}
}
bool WellInterfaceGeneric::isOperableAndSolvable() const
template<class Scalar>
bool WellInterfaceGeneric<Scalar>::isOperableAndSolvable() const
{
return operability_status_.isOperableAndSolvable();
}
bool WellInterfaceGeneric::useVfpExplicit() const
template<class Scalar>
bool WellInterfaceGeneric<Scalar>::useVfpExplicit() const
{
const auto& wvfpexp = well_ecl_.getWVFPEXP();
return (wvfpexp.explicit_lookup() || operability_status_.use_vfpexplicit);
}
bool WellInterfaceGeneric::thpLimitViolatedButNotSwitched() const
template<class Scalar>
bool WellInterfaceGeneric<Scalar>::thpLimitViolatedButNotSwitched() const
{
return operability_status_.thp_limit_violated_but_not_switched;
}
double WellInterfaceGeneric::getALQ(const WellState<double>& well_state) const
template<class Scalar>
Scalar WellInterfaceGeneric<Scalar>::
getALQ(const WellState<Scalar>& well_state) const
{
// no alq for injectors.
if (isInjector())
@ -514,8 +565,10 @@ double WellInterfaceGeneric::getALQ(const WellState<double>& well_state) const
return well_state.getALQ(name());
}
void WellInterfaceGeneric::reportWellSwitching(const SingleWellState<double> &ws,
DeferredLogger& deferred_logger) const
template<class Scalar>
void WellInterfaceGeneric<Scalar>::
reportWellSwitching(const SingleWellState<Scalar> &ws,
DeferredLogger& deferred_logger) const
{
if (well_control_log_.empty())
return;
@ -534,7 +587,9 @@ void WellInterfaceGeneric::reportWellSwitching(const SingleWellState<double> &ws
}
}
bool WellInterfaceGeneric::isPressureControlled(const WellState<double>& well_state) const
template<class Scalar>
bool WellInterfaceGeneric<Scalar>::
isPressureControlled(const WellState<Scalar>& well_state) const
{
const auto& ws = well_state.well(this->index_of_well_);
if (this->isInjector()) {
@ -548,10 +603,10 @@ bool WellInterfaceGeneric::isPressureControlled(const WellState<double>& well_st
}
}
bool WellInterfaceGeneric::wellUnderZeroRateTarget(const SummaryState& summary_state,
const WellState<double>& well_state) const
template<class Scalar>
bool WellInterfaceGeneric<Scalar>::
wellUnderZeroRateTarget(const SummaryState& summary_state,
const WellState<Scalar>& well_state) const
{
if (this->isProducer()) { // producers
const auto prod_controls = this->well_ecl_.productionControls(summary_state);
@ -564,25 +619,29 @@ bool WellInterfaceGeneric::wellUnderZeroRateTarget(const SummaryState& summary_s
}
}
bool WellInterfaceGeneric::stopppedOrZeroRateTarget(const SummaryState& summary_state,
const WellState<double>& well_state) const
template<class Scalar>
bool WellInterfaceGeneric<Scalar>::
stopppedOrZeroRateTarget(const SummaryState& summary_state,
const WellState<Scalar>& well_state) const
{
return (this->wellIsStopped() || this->wellUnderZeroRateTarget(summary_state, well_state));
}
void WellInterfaceGeneric::resetWellOperability()
template<class Scalar>
void WellInterfaceGeneric<Scalar>::resetWellOperability()
{
this->operability_status_.resetOperability();
}
double WellInterfaceGeneric::wmicrobes_() const
template<class Scalar>
Scalar WellInterfaceGeneric<Scalar>::wmicrobes_() const
{
auto injectorType = this->well_ecl_.injectorType();
if (injectorType == InjectorType::WATER) {
WellMICPProperties microbes = this->well_ecl_.getMICPProperties();
const double microbial_injection_concentration = microbes.m_microbialConcentration;
const Scalar microbial_injection_concentration = microbes.m_microbialConcentration;
return microbial_injection_concentration;
} else {
// Not a water injection well => no microbes.
@ -590,7 +649,8 @@ double WellInterfaceGeneric::wmicrobes_() const
}
}
double WellInterfaceGeneric::wfoam_() const
template<class Scalar>
Scalar WellInterfaceGeneric<Scalar>::wfoam_() const
{
auto injectorType = this->well_ecl_.injectorType();
@ -603,7 +663,8 @@ double WellInterfaceGeneric::wfoam_() const
}
}
double WellInterfaceGeneric::wsalt_() const
template<class Scalar>
Scalar WellInterfaceGeneric<Scalar>::wsalt_() const
{
auto injectorType = this->well_ecl_.injectorType();
@ -616,13 +677,14 @@ double WellInterfaceGeneric::wsalt_() const
}
}
double WellInterfaceGeneric::woxygen_() const
template<class Scalar>
Scalar WellInterfaceGeneric<Scalar>::woxygen_() const
{
auto injectorType = this->well_ecl_.injectorType();
if (injectorType == InjectorType::WATER) {
WellMICPProperties oxygen = this->well_ecl_.getMICPProperties();
const double oxygen_injection_concentration = oxygen.m_oxygenConcentration;
const Scalar oxygen_injection_concentration = oxygen.m_oxygenConcentration;
return oxygen_injection_concentration;
} else {
// Not a water injection well => no oxygen.
@ -630,13 +692,14 @@ double WellInterfaceGeneric::woxygen_() const
}
}
double WellInterfaceGeneric::wpolymer_() const
template<class Scalar>
Scalar WellInterfaceGeneric<Scalar>::wpolymer_() const
{
auto injectorType = this->well_ecl_.injectorType();
if (injectorType == InjectorType::WATER) {
WellPolymerProperties polymer = this->well_ecl_.getPolymerProperties();
const double polymer_injection_concentration = polymer.m_polymerConcentration;
const Scalar polymer_injection_concentration = polymer.m_polymerConcentration;
return polymer_injection_concentration;
} else {
// Not a water injection well => no polymer.
@ -644,13 +707,14 @@ double WellInterfaceGeneric::wpolymer_() const
}
}
double WellInterfaceGeneric::wurea_() const
template<class Scalar>
Scalar WellInterfaceGeneric<Scalar>::wurea_() const
{
auto injectorType = this->well_ecl_.injectorType();
if (injectorType == InjectorType::WATER) {
WellMICPProperties urea = this->well_ecl_.getMICPProperties();
const double urea_injection_concentration = urea.m_ureaConcentration / 10.; //Dividing by scaling factor 10
const Scalar urea_injection_concentration = urea.m_ureaConcentration / 10.; //Dividing by scaling factor 10
return urea_injection_concentration;
} else {
// Not a water injection well => no urea.
@ -658,24 +722,28 @@ double WellInterfaceGeneric::wurea_() const
}
}
int WellInterfaceGeneric::polymerTable_() const
template<class Scalar>
int WellInterfaceGeneric<Scalar>::polymerTable_() const
{
return this->well_ecl_.getPolymerProperties().m_skprpolytable;
}
int WellInterfaceGeneric::polymerWaterTable_() const
template<class Scalar>
int WellInterfaceGeneric<Scalar>::polymerWaterTable_() const
{
return this->well_ecl_.getPolymerProperties().m_skprwattable;
}
int WellInterfaceGeneric::polymerInjTable_() const
template<class Scalar>
int WellInterfaceGeneric<Scalar>::polymerInjTable_() const
{
return this->well_ecl_.getPolymerProperties().m_plymwinjtable;
}
std::pair<bool,bool> WellInterfaceGeneric::
computeWellPotentials(std::vector<double>& well_potentials,
const WellState<double>& well_state)
template<class Scalar>
std::pair<bool,bool> WellInterfaceGeneric<Scalar>::
computeWellPotentials(std::vector<Scalar>& well_potentials,
const WellState<Scalar>& well_state)
{
const int np = this->number_of_phases_;
well_potentials.resize(np, 0.0);
@ -711,8 +779,8 @@ computeWellPotentials(std::vector<double>& well_potentials,
if (!this->changed_to_open_this_step_ &&
(thp_controlled_well || bhp_controlled_well)) {
double total_rate = 0.0;
const double sign = this->isInjector() ? 1.0 : -1.0;
Scalar total_rate = 0.0;
const Scalar sign = this->isInjector() ? 1.0 : -1.0;
for (int phase = 0; phase < np; ++phase){
total_rate += sign * ws.surface_rates[phase];
}
@ -730,13 +798,14 @@ computeWellPotentials(std::vector<double>& well_potentials,
return {compute_potential, bhp_controlled_well};
}
void WellInterfaceGeneric::
checkNegativeWellPotentials(std::vector<double>& well_potentials,
template<class Scalar>
void WellInterfaceGeneric<Scalar>::
checkNegativeWellPotentials(std::vector<Scalar>& well_potentials,
const bool checkOperability,
DeferredLogger& deferred_logger)
{
const double sign = this->isInjector() ? 1.0 : -1.0;
double total_potential = 0.0;
const Scalar sign = this->isInjector() ? 1.0 : -1.0;
Scalar total_potential = 0.0;
for (int phase = 0; phase < this->number_of_phases_; ++phase) {
well_potentials[phase] *= sign;
total_potential += well_potentials[phase];
@ -750,9 +819,10 @@ checkNegativeWellPotentials(std::vector<double>& well_potentials,
}
}
void WellInterfaceGeneric::
template<class Scalar>
void WellInterfaceGeneric<Scalar>::
prepareForPotentialCalculations(const SummaryState& summary_state,
WellState<double>& well_state,
WellState<Scalar>& well_state,
Well::InjectionControls& inj_controls,
Well::ProductionControls& prod_controls) const
{
@ -793,4 +863,6 @@ prepareForPotentialCalculations(const SummaryState& summary_state,
}
}
template class WellInterfaceGeneric<double>;
} // namespace Opm

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

@ -47,6 +47,7 @@ template<class Scalar> class SingleWellState;
class Group;
class Schedule;
template<class Scalar>
class WellInterfaceGeneric {
public:
WellInterfaceGeneric(const Well& well,
@ -76,7 +77,7 @@ public:
/// Index of well in the wells struct and wellState
int indexOfWell() const;
void adaptRatesForVFP(std::vector<double>& rates) const;
void adaptRatesForVFP(std::vector<Scalar>& rates) const;
const Well& wellEcl() const;
Well& wellEcl();
@ -94,127 +95,91 @@ public:
void closeCompletions(const WellTestState& wellTestState);
void setVFPProperties(const VFPProperties* vfp_properties_arg);
void setPrevSurfaceRates(WellState<double>& well_state,
const WellState<double>& prev_well_state) const;
void setPrevSurfaceRates(WellState<Scalar>& well_state,
const WellState<Scalar>& prev_well_state) const;
void setGuideRate(const GuideRate* guide_rate_arg);
void setWellEfficiencyFactor(const double efficiency_factor);
void setWellEfficiencyFactor(const Scalar efficiency_factor);
void setRepRadiusPerfLength();
void setWsolvent(const double wsolvent);
void setDynamicThpLimit(const double thp_limit);
std::optional<double> getDynamicThpLimit() const;
void setWsolvent(const Scalar wsolvent);
void setDynamicThpLimit(const Scalar thp_limit);
std::optional<Scalar> getDynamicThpLimit() const;
void updatePerforatedCell(std::vector<bool>& is_cell_perforated);
/// Returns true if the well has one or more THP limits/constraints.
bool wellHasTHPConstraints(const SummaryState& summaryState) const;
void stopWell() {
this->wellStatus_ = Well::Status::STOP;
}
void stopWell() { this->wellStatus_ = Well::Status::STOP; }
void openWell() { this->wellStatus_ = Well::Status::OPEN; }
void openWell() {
this->wellStatus_ = Well::Status::OPEN;
}
bool wellIsStopped() const { return this->wellStatus_ == Well::Status::STOP; }
bool wellIsStopped() const {
return this->wellStatus_ == Well::Status::STOP;
}
int currentStep() const { return this->current_step_; }
int currentStep() const {
return this->current_step_;
}
int pvtRegionIdx() const { return pvtRegionIdx_; }
int pvtRegionIdx() const {
return pvtRegionIdx_;
}
const GuideRate* guideRate() const { return guide_rate_; }
const GuideRate* guideRate() const {
return guide_rate_;
}
int numComponents() const { return num_components_; }
int numComponents() const {
return num_components_;
}
int numPhases() const { return number_of_phases_; }
int numPhases() const {
return number_of_phases_;
}
int numPerfs() const { return number_of_perforations_; }
int numPerfs() const {
return number_of_perforations_;
}
Scalar refDepth() const { return ref_depth_; }
double refDepth() const {
return ref_depth_;
}
Scalar gravity() const { return gravity_; }
double gravity() const {
return gravity_;
}
const VFPProperties* vfpProperties() const { return vfp_properties_; }
const VFPProperties* vfpProperties() const {
return vfp_properties_;
}
const ParallelWellInfo& parallelWellInfo() const { return parallel_well_info_; }
const ParallelWellInfo& parallelWellInfo() const {
return parallel_well_info_;
}
const std::vector<Scalar>& perfDepth() const { return perf_depth_; }
const std::vector<double>& perfDepth() const {
return perf_depth_;
}
std::vector<Scalar>& perfDepth() { return perf_depth_; }
std::vector<double>& perfDepth() {
return perf_depth_;
}
const std::vector<Scalar>& wellIndex() const { return well_index_; }
const std::vector<double>& wellIndex() const {
return well_index_;
}
const std::map<int,std::vector<int>>& getCompletions() const { return completions_; }
const std::map<int,std::vector<int>>& getCompletions() const {
return completions_;
}
double getTHPConstraint(const SummaryState& summaryState) const;
double getALQ(const WellState<double>& well_state) const;
double wsolvent() const;
double rsRvInj() const;
Scalar getTHPConstraint(const SummaryState& summaryState) const;
Scalar getALQ(const WellState<Scalar>& well_state) const;
Scalar wsolvent() const;
Scalar rsRvInj() const;
// at the beginning of the time step, we check what inj_multiplier from the previous running
void initInjMult(const std::vector<double>& max_inj_mult);
void initInjMult(const std::vector<Scalar>& max_inj_mult);
// update the InjMult information at the end of the time step, so it can be used for later.
void updateInjMult(std::vector<double>& inj_multipliers, DeferredLogger& deferred_logger) const;
void updateInjMult(std::vector<Scalar>& inj_multipliers,
DeferredLogger& deferred_logger) const;
// Note:: for multisegment wells, bhp is actually segment pressure in practice based on observation
// it might change in the future
double getInjMult(const int perf, const double bhp, const double perf_pres) const;
Scalar getInjMult(const int perf, const Scalar bhp, const Scalar perf_pres) const;
// whether a well is specified with a non-zero and valid VFP table number
bool isVFPActive(DeferredLogger& deferred_logger) const;
void reportWellSwitching(const SingleWellState<double>& ws, DeferredLogger& deferred_logger) const;
void reportWellSwitching(const SingleWellState<Scalar>& ws,
DeferredLogger& deferred_logger) const;
bool changedToOpenThisStep() const {
return this->changed_to_open_this_step_;
}
bool changedToOpenThisStep() const { return this->changed_to_open_this_step_; }
void updateWellTestState(const SingleWellState<double>& ws,
void updateWellTestState(const SingleWellState<Scalar>& ws,
const double& simulationTime,
const bool& writeMessageToOPMLog,
WellTestState& wellTestState,
DeferredLogger& deferred_logger) const;
bool isPressureControlled(const WellState<double>& well_state) const;
bool isPressureControlled(const WellState<Scalar>& well_state) const;
bool stopppedOrZeroRateTarget(const SummaryState& summary_state,
const WellState<double>& well_state) const;
const WellState<Scalar>& well_state) const;
double wellEfficiencyFactor() const
{ return well_efficiency_factor_; }
Scalar wellEfficiencyFactor() const { return well_efficiency_factor_; }
//! \brief Update filter cake multipliers.
void updateFilterCakeMultipliers(const std::vector<double>& inj_fc_multiplier)
void updateFilterCakeMultipliers(const std::vector<Scalar>& inj_fc_multiplier)
{
inj_fc_multiplier_ = inj_fc_multiplier;
}
@ -224,36 +189,38 @@ public:
protected:
bool getAllowCrossFlow() const;
double wmicrobes_() const;
double wfoam_() const;
double woxygen_() const;
double wpolymer_() const;
double wsalt_() const;
double wurea_() const;
Scalar wmicrobes_() const;
Scalar wfoam_() const;
Scalar woxygen_() const;
Scalar wpolymer_() const;
Scalar wsalt_() const;
Scalar wurea_() const;
int polymerTable_() const;
int polymerInjTable_() const;
int polymerWaterTable_() const;
bool wellUnderZeroRateTarget(const SummaryState& summary_state,
const WellState<double>& well_state) const;
const WellState<Scalar>& well_state) const;
std::pair<bool,bool>
computeWellPotentials(std::vector<double>& well_potentials,
const WellState<double>& well_state);
computeWellPotentials(std::vector<Scalar>& well_potentials,
const WellState<Scalar>& well_state);
void checkNegativeWellPotentials(std::vector<double>& well_potentials,
void checkNegativeWellPotentials(std::vector<Scalar>& well_potentials,
const bool checkOperability,
DeferredLogger& deferred_logger);
void prepareForPotentialCalculations(const SummaryState& summary_state,
WellState<double>& well_state,
WellState<Scalar>& well_state,
Well::InjectionControls& inj_controls,
Well::ProductionControls& prod_controls) const;
// definition of the struct OperabilityStatus
struct OperabilityStatus {
bool isOperableAndSolvable() const {
struct OperabilityStatus
{
bool isOperableAndSolvable() const
{
if (!operable_under_only_bhp_limit || !solvable || has_negative_potentials) {
return false;
} else {
@ -261,15 +228,18 @@ protected:
}
}
bool isOperableUnderBHPLimit() const {
bool isOperableUnderBHPLimit() const
{
return operable_under_only_bhp_limit && obey_thp_limit_under_bhp_limit;
}
bool isOperableUnderTHPLimit() const {
bool isOperableUnderTHPLimit() const
{
return can_obtain_bhp_with_thp_limit && obey_bhp_limit_with_thp_limit;
}
void resetOperability() {
void resetOperability()
{
operable_under_only_bhp_limit = true;
obey_thp_limit_under_bhp_limit = true;
can_obtain_bhp_with_thp_limit = true;
@ -322,29 +292,29 @@ protected:
// Q = IPR_A - BHP * IPR_B
// TODO: it minght need to go to WellInterface, let us implement it in StandardWell first
// it is only updated and used for producers for now
mutable std::vector<double> ipr_a_;
mutable std::vector<double> ipr_b_;
mutable std::vector<Scalar> ipr_a_;
mutable std::vector<Scalar> ipr_b_;
// cell index for each well perforation
std::vector<int> well_cells_;
// well index for each perforation
std::vector<double> well_index_;
std::vector<Scalar> well_index_;
// number of the perforations for this well
int number_of_perforations_;
// depth for each perforation
std::vector<double> perf_depth_;
std::vector<Scalar> perf_depth_;
// representative radius of the perforations, used in shear calculation
std::vector<double> perf_rep_radius_;
std::vector<Scalar> perf_rep_radius_;
// length of the perforations, use in shear calculation
std::vector<double> perf_length_;
std::vector<Scalar> perf_length_;
// well bore diameter
std::vector<double> bore_diameters_;
std::vector<Scalar> bore_diameters_;
/*
* completions_ contains the mapping from completion id to connection indices
@ -364,7 +334,7 @@ protected:
std::map<int, std::vector<int>> completions_;
// reference depth for the BHP
double ref_depth_;
Scalar ref_depth_;
// saturation table nubmer for each well perforation
std::vector<int> saturation_table_number_;
@ -373,25 +343,25 @@ protected:
const PhaseUsage* phase_usage_;
double gravity_;
double wsolvent_;
std::optional<double> dynamic_thp_limit_;
Scalar gravity_;
Scalar wsolvent_;
std::optional<Scalar> dynamic_thp_limit_;
// recording the multiplier calculate from the keyword WINJMULT during the time step
mutable std::vector<double> inj_multiplier_;
mutable std::vector<Scalar> inj_multiplier_;
// the injection multiplier from the previous running, it is mostly used for CIRR mode
// which intends to keep the fracturing open
std::vector<double> prev_inj_multiplier_;
std::vector<Scalar> prev_inj_multiplier_;
// the multiplier due to injection filtration cake
std::vector<double> inj_fc_multiplier_;
std::vector<Scalar> inj_fc_multiplier_;
double well_efficiency_factor_;
Scalar well_efficiency_factor_;
const VFPProperties* vfp_properties_;
const GuideRate* guide_rate_;
std::vector< std::string> well_control_log_;
std::vector<std::string> well_control_log_;
bool changed_to_open_this_step_ = true;
};

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

@ -91,7 +91,7 @@ modelCompIdxToFlowCompIdx(const unsigned compIdx) const
}
template<class FluidSystem, class Indices>
double
typename WellInterfaceIndices<FluidSystem,Indices>::Scalar
WellInterfaceIndices<FluidSystem,Indices>::
scalingFactor(const int phaseIdx) const
{

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

@ -27,8 +27,7 @@
#include <opm/simulators/wells/WellInterfaceFluidSystem.hpp>
namespace Opm
{
namespace Opm {
template<class FluidSystem, class Indices>
class WellInterfaceIndices : public WellInterfaceFluidSystem<FluidSystem>
@ -42,7 +41,7 @@ public:
int flowPhaseToModelCompIdx(const int phaseIdx) const;
int modelCompIdxToFlowCompIdx(const unsigned compIdx) const;
double scalingFactor(const int phaseIdx) const;
Scalar scalingFactor(const int phaseIdx) const;
template <class EvalWell>
Eval restrictEval(const EvalWell& in) const

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

@ -89,10 +89,10 @@ namespace Opm
void
WellInterface<TypeTag>::
init(const PhaseUsage* phase_usage_arg,
const std::vector<double>& /* depth_arg */,
const double gravity_arg,
const std::vector<Scalar>& /* depth_arg */,
const Scalar 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;
@ -105,7 +105,7 @@ namespace Opm
template<typename TypeTag>
double
typename WellInterface<TypeTag>::Scalar
WellInterface<TypeTag>::
wpolymer() const
{
@ -121,7 +121,7 @@ namespace Opm
template<typename TypeTag>
double
typename WellInterface<TypeTag>::Scalar
WellInterface<TypeTag>::
wfoam() const
{
@ -135,7 +135,7 @@ namespace Opm
template<typename TypeTag>
double
typename WellInterface<TypeTag>::Scalar
WellInterface<TypeTag>::
wsalt() const
{
@ -147,7 +147,7 @@ namespace Opm
}
template<typename TypeTag>
double
typename WellInterface<TypeTag>::Scalar
WellInterface<TypeTag>::
wmicrobes() const
{
@ -159,7 +159,7 @@ namespace Opm
}
template<typename TypeTag>
double
typename WellInterface<TypeTag>::Scalar
WellInterface<TypeTag>::
woxygen() const
{
@ -177,7 +177,7 @@ namespace Opm
// vtkblackoilmicpmodule.hh and ecloutputblackoilmodel.hh respectively).
template<typename TypeTag>
double
typename WellInterface<TypeTag>::Scalar
WellInterface<TypeTag>::
wurea() const
{
@ -273,7 +273,7 @@ namespace Opm
const GroupState<Scalar>& group_state,
const Well::InjectionControls& inj_controls,
const Well::ProductionControls& prod_controls,
const double wqTotal,
const Scalar wqTotal,
DeferredLogger& deferred_logger,
const bool fixed_control,
const bool fixed_status)
@ -285,7 +285,7 @@ namespace Opm
return false;
}
const double sgn = this->isInjector() ? 1.0 : -1.0;
const Scalar sgn = this->isInjector() ? 1.0 : -1.0;
if (!this->wellIsStopped()){
if (wqTotal*sgn <= 0.0 && !fixed_status){
this->stopWell();
@ -295,7 +295,7 @@ namespace Opm
if (!fixed_control) {
auto& ws = well_state.well(this->index_of_well_);
const bool hasGroupControl = this->isInjector() ? inj_controls.hasControl(Well::InjectorCMode::GRUP) :
prod_controls.hasControl(Well::ProducerCMode::GRUP);
prod_controls.hasControl(Well::ProducerCMode::GRUP);
changed = this->checkIndividualConstraints(ws, summary_state, deferred_logger, inj_controls, prod_controls);
if (hasGroupControl) {
@ -318,22 +318,32 @@ namespace Opm
}
} else if (!fixed_status){
// well is stopped, check if current bhp allows reopening
const double bhp = well_state.well(this->index_of_well_).bhp;
double prod_limit = prod_controls.bhp_limit;
double inj_limit = inj_controls.bhp_limit;
const Scalar bhp = well_state.well(this->index_of_well_).bhp;
Scalar prod_limit = prod_controls.bhp_limit;
Scalar inj_limit = inj_controls.bhp_limit;
const bool has_thp = this->wellHasTHPConstraints(summary_state);
if (has_thp){
std::vector<double> rates(this->num_components_);
std::vector<Scalar> rates(this->num_components_);
if (this->isInjector()){
const double bhp_thp = WellBhpThpCalculator(*this).calculateBhpFromThp(well_state, rates, this->well_ecl_, summary_state, this->getRefDensity(), deferred_logger);
inj_limit = std::min(bhp_thp, inj_controls.bhp_limit);
const Scalar bhp_thp = WellBhpThpCalculator(*this).
calculateBhpFromThp(well_state, rates,
this->well_ecl_,
summary_state,
this->getRefDensity(),
deferred_logger);
inj_limit = std::min(bhp_thp, static_cast<Scalar>(inj_controls.bhp_limit));
} else {
// if the well can operate, it must at least be able to produce at the lowest bhp of the bhp-curve (explicit fractions)
const double bhp_min = WellBhpThpCalculator(*this).calculateMinimumBhpFromThp(well_state, this->well_ecl_, summary_state, this->getRefDensity());
prod_limit = std::max(bhp_min, prod_controls.bhp_limit);
// if the well can operate, it must at least be able to produce
// at the lowest bhp of the bhp-curve (explicit fractions)
const Scalar bhp_min = WellBhpThpCalculator(*this).
calculateMinimumBhpFromThp(well_state,
this->well_ecl_,
summary_state,
this->getRefDensity());
prod_limit = std::max(bhp_min, static_cast<Scalar>(prod_controls.bhp_limit));
}
}
const double bhp_diff = (this->isInjector())? inj_limit - bhp: bhp - prod_limit;
const Scalar bhp_diff = (this->isInjector())? inj_limit - bhp: bhp - prod_limit;
if (bhp_diff > 0){
this->openWell();
well_state.well(this->index_of_well_).bhp = (this->isInjector())? inj_limit : prod_limit;
@ -401,19 +411,25 @@ namespace Opm
deferred_logger.debug(msg);
return;
}
std::vector<double> potentials;
std::vector<Scalar> potentials;
try {
computeWellPotentials(simulator, well_state_copy, potentials, deferred_logger);
} catch (const std::exception& e) {
const std::string msg = std::string("well ") + this->name() + std::string(": computeWellPotentials() failed during testing for re-opening: ") + e.what();
const std::string msg = fmt::format("well {}: computeWellPotentials() "
"failed during testing for re-opening: ",
this->name(), e.what());
deferred_logger.info(msg);
return;
}
const int np = well_state_copy.numPhases();
for (int p = 0; p < np; ++p) {
ws.well_potentials[p] = std::max(0.0, potentials[p]);
ws.well_potentials[p] = std::max(Scalar{0.0}, potentials[p]);
}
this->updateWellTestState(well_state_copy.well(this->indexOfWell()), simulation_time, /*writeMessageToOPMLog=*/ false, welltest_state_temp, deferred_logger);
this->updateWellTestState(well_state_copy.well(this->indexOfWell()),
simulation_time,
/*writeMessageToOPMLog=*/ false,
welltest_state_temp,
deferred_logger);
this->closeCompletions(welltest_state_temp);
// Stop testing if the well is closed or shut due to all completions shut
@ -510,7 +526,8 @@ namespace Opm
} else {
// solve well with the estimated target bhp (or limit)
ws.thp = this->getTHPConstraint(summary_state);
const double bhp = std::max(bhp_target.value(), prod_controls.bhp_limit);
const Scalar bhp = std::max(bhp_target.value(),
static_cast<Scalar>(prod_controls.bhp_limit));
solveWellWithBhp(simulator, dt, bhp, well_state, deferred_logger);
}
}
@ -531,8 +548,8 @@ namespace Opm
this->operability_status_.use_vfpexplicit = true;
auto bhp_stable = WellBhpThpCalculator(*this).estimateStableBhp(well_state, this->well_ecl_, rates, this->getRefDensity(), summary_state);
// if we find an intersection with a sufficiently lower bhp, re-solve equations
const double reltol = 1e-3;
const double cur_bhp = ws.bhp;
const Scalar reltol = 1e-3;
const Scalar cur_bhp = ws.bhp;
if (bhp_stable.has_value() && cur_bhp - bhp_stable.value() > cur_bhp*reltol){
const auto msg = fmt::format("Well {} converged to an unstable solution, re-solving", this->name());
deferred_logger.debug(msg);
@ -557,10 +574,16 @@ namespace Opm
this->stopWell();
} else {
// solve well with the estimated target bhp (or limit)
const double bhp = std::max(bhp_target.value(), prod_controls.bhp_limit);
const Scalar bhp = std::max(bhp_target.value(),
static_cast<Scalar>(prod_controls.bhp_limit));
solveWellWithBhp(simulator, dt, bhp, well_state, deferred_logger);
ws.thp = this->getTHPConstraint(summary_state);
converged = this->iterateWellEqWithSwitching(simulator, dt, inj_controls, prod_controls, well_state, group_state, deferred_logger);
converged = this->iterateWellEqWithSwitching(simulator, dt,
inj_controls,
prod_controls,
well_state,
group_state,
deferred_logger);
}
}
// update operability
@ -571,7 +594,7 @@ namespace Opm
}
template<typename TypeTag>
std::optional<double>
std::optional<typename WellInterface<TypeTag>::Scalar>
WellInterface<TypeTag>::
estimateOperableBhp(const Simulator& simulator,
const double dt,
@ -581,7 +604,7 @@ namespace Opm
{
// Given an unconverged well or closed well, estimate an operable bhp (if any)
// Get minimal bhp from vfp-curve
double bhp_min = WellBhpThpCalculator(*this).calculateMinimumBhpFromThp(well_state, this->well_ecl_, summary_state, this->getRefDensity());
Scalar bhp_min = WellBhpThpCalculator(*this).calculateMinimumBhpFromThp(well_state, this->well_ecl_, summary_state, this->getRefDensity());
// Solve
const bool converged = solveWellWithBhp(simulator, dt, bhp_min, well_state, deferred_logger);
if (!converged || this->wellIsStopped()) {
@ -598,7 +621,7 @@ namespace Opm
WellInterface<TypeTag>::
solveWellWithBhp(const Simulator& simulator,
const double dt,
const double bhp,
const Scalar bhp,
WellState<Scalar>& well_state,
DeferredLogger& deferred_logger)
{
@ -750,9 +773,7 @@ namespace Opm
const GroupState<Scalar>& group_state,
DeferredLogger& deferred_logger)
{
prepareWellBeforeAssembling(simulator, dt, well_state, group_state, deferred_logger);
assembleWellEqWithoutIteration(simulator, dt, well_state, group_state, deferred_logger);
}
@ -805,18 +826,20 @@ namespace Opm
}
if (this->operability_status_.has_negative_potentials) {
auto well_state_copy = well_state;
std::vector<double> potentials;
std::vector<Scalar> potentials;
try {
computeWellPotentials(simulator, well_state_copy, potentials, deferred_logger);
} catch (const std::exception& e) {
const std::string msg = std::string("well ") + this->name() + std::string(": computeWellPotentials() failed during attempt to recompute potentials for well : ") + e.what();
const std::string msg = fmt::format("well {}: computeWellPotentials() failed "
"during attempt to recompute potentials for well: ",
this->name(), e.what());
deferred_logger.info(msg);
this->operability_status_.has_negative_potentials = true;
}
auto& ws = well_state.well(this->indexOfWell());
const int np = well_state.numPhases();
for (int p = 0; p < np; ++p) {
ws.well_potentials[p] = std::max(0.0, potentials[p]);
ws.well_potentials[p] = std::max(Scalar{0.0}, potentials[p]);
}
}
this->changed_to_open_this_step_ = false;
@ -852,7 +875,8 @@ namespace Opm
template<typename TypeTag>
typename WellInterface<TypeTag>::Scalar
WellInterface<TypeTag>::volumetricSurfaceRateForConnection(int cellIdx, int phaseIdx) const {
WellInterface<TypeTag>::volumetricSurfaceRateForConnection(int cellIdx, int phaseIdx) const
{
for (int perfIdx = 0; perfIdx < this->number_of_perforations_; ++perfIdx) {
if (this->cells()[perfIdx] == cellIdx) {
const unsigned activeCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx));
@ -896,12 +920,11 @@ namespace Opm
template<typename TypeTag>
bool
WellInterface<TypeTag>::
gliftBeginTimeStepWellTestIterateWellEquations(
const Simulator& simulator,
const double dt,
WellState<Scalar>& well_state,
const GroupState<Scalar>& group_state,
DeferredLogger& deferred_logger)
gliftBeginTimeStepWellTestIterateWellEquations(const Simulator& simulator,
const double dt,
WellState<Scalar>& well_state,
const GroupState<Scalar>& group_state,
DeferredLogger& deferred_logger)
{
const auto& well_name = this->name();
assert(this->wellHasTHPConstraints(simulator.vanguard().summaryState()));
@ -954,7 +977,7 @@ namespace Opm
}
const auto& gl_well = glo.well(well_name);
auto& max_alq_optional = gl_well.max_rate();
double max_alq;
Scalar max_alq;
if (max_alq_optional) {
max_alq = *max_alq_optional;
}
@ -1080,7 +1103,7 @@ namespace Opm
const auto current = ws.injection_cmode;
switch(current) {
switch (current) {
case Well::InjectorCMode::RATE:
{
ws.surface_rates[phasePos] = (1.0 - this->rsRvInj()) * controls.surface_rate;
@ -1098,9 +1121,9 @@ namespace Opm
case Well::InjectorCMode::RESV:
{
std::vector<double> convert_coeff(this->number_of_phases_, 1.0);
std::vector<Scalar> convert_coeff(this->number_of_phases_, 1.0);
this->rateConverter_.calcCoeff(/*fipreg*/ 0, this->pvtRegionIdx_, convert_coeff);
const double coeff = convert_coeff[phasePos];
const Scalar coeff = convert_coeff[phasePos];
ws.surface_rates[phasePos] = controls.reservoir_rate/coeff;
break;
}
@ -1108,7 +1131,7 @@ namespace Opm
case Well::InjectorCMode::THP:
{
auto rates = ws.surface_rates;
double bhp = WellBhpThpCalculator(*this).calculateBhpFromThp(well_state,
Scalar bhp = WellBhpThpCalculator(*this).calculateBhpFromThp(well_state,
rates,
well,
summaryState,
@ -1120,7 +1143,7 @@ namespace Opm
// if the total rates are negative or zero
// we try to provide a better intial well rate
// using the well potentials
double total_rate = std::accumulate(rates.begin(), rates.end(), 0.0);
Scalar total_rate = std::accumulate(rates.begin(), rates.end(), 0.0);
if (total_rate <= 0.0)
ws.surface_rates = ws.well_potentials;
@ -1129,7 +1152,7 @@ namespace Opm
case Well::InjectorCMode::BHP:
{
ws.bhp = controls.bhp_limit;
double total_rate = 0.0;
Scalar total_rate = 0.0;
for (int p = 0; p<np; ++p) {
total_rate += ws.surface_rates[p];
}
@ -1145,8 +1168,8 @@ namespace Opm
{
assert(well.isAvailableForGroupControl());
const auto& group = schedule.getGroup(well.groupName(), this->currentStep());
const double efficiencyFactor = well.getEfficiencyFactor();
std::optional<double> target =
const Scalar efficiencyFactor = well.getEfficiencyFactor();
std::optional<Scalar> target =
this->getGroupInjectionTargetRate(group,
well_state,
group_state,
@ -1171,7 +1194,7 @@ namespace Opm
// the zero rate target, then we can use to retain the composition information
// within the wellbore from the previous result, and hopefully it is a good
// initial guess for the zero rate target.
ws.surface_rates[phasePos] = std::max(1.e-7, ws.surface_rates[phasePos]);
ws.surface_rates[phasePos] = std::max(Scalar{1.e-7}, ws.surface_rates[phasePos]);
if (ws.bhp == 0.) {
ws.bhp = controls.bhp_limit;
@ -1185,7 +1208,7 @@ namespace Opm
switch (current) {
case Well::ProducerCMode::ORAT:
{
double current_rate = -ws.surface_rates[ pu.phase_pos[Oil] ];
Scalar current_rate = -ws.surface_rates[ pu.phase_pos[Oil] ];
// for trivial rates or opposite direction we don't just scale the rates
// but use either the potentials or the mobility ratio to initial the well rates
if (current_rate > 0.0) {
@ -1193,7 +1216,7 @@ namespace Opm
ws.surface_rates[p] *= controls.oil_rate/current_rate;
}
} else {
const std::vector<double> fractions = initialWellRateFractions(simulator, well_state);
const std::vector<Scalar> fractions = initialWellRateFractions(simulator, well_state);
double control_fraction = fractions[pu.phase_pos[Oil]];
if (control_fraction != 0.0) {
for (int p = 0; p<np; ++p) {
@ -1205,7 +1228,7 @@ namespace Opm
}
case Well::ProducerCMode::WRAT:
{
double current_rate = -ws.surface_rates[ pu.phase_pos[Water] ];
Scalar current_rate = -ws.surface_rates[ pu.phase_pos[Water] ];
// for trivial rates or opposite direction we don't just scale the rates
// but use either the potentials or the mobility ratio to initial the well rates
if (current_rate > 0.0) {
@ -1213,11 +1236,11 @@ namespace Opm
ws.surface_rates[p] *= controls.water_rate/current_rate;
}
} else {
const std::vector<double> fractions = initialWellRateFractions(simulator, well_state);
double control_fraction = fractions[pu.phase_pos[Water]];
const std::vector<Scalar> fractions = initialWellRateFractions(simulator, well_state);
const Scalar control_fraction = fractions[pu.phase_pos[Water]];
if (control_fraction != 0.0) {
for (int p = 0; p<np; ++p) {
ws.surface_rates[p] = - fractions[p] * controls.water_rate/control_fraction;
ws.surface_rates[p] = - fractions[p] * controls.water_rate / control_fraction;
}
}
}
@ -1225,7 +1248,7 @@ namespace Opm
}
case Well::ProducerCMode::GRAT:
{
double current_rate = -ws.surface_rates[pu.phase_pos[Gas] ];
Scalar current_rate = -ws.surface_rates[pu.phase_pos[Gas] ];
// or trivial rates or opposite direction we don't just scale the rates
// but use either the potentials or the mobility ratio to initial the well rates
if (current_rate > 0.0) {
@ -1233,11 +1256,11 @@ namespace Opm
ws.surface_rates[p] *= controls.gas_rate/current_rate;
}
} else {
const std::vector<double> fractions = initialWellRateFractions(simulator, well_state);
double control_fraction = fractions[pu.phase_pos[Gas]];
const std::vector<Scalar > fractions = initialWellRateFractions(simulator, well_state);
const Scalar control_fraction = fractions[pu.phase_pos[Gas]];
if (control_fraction != 0.0) {
for (int p = 0; p<np; ++p) {
ws.surface_rates[p] = - fractions[p] * controls.gas_rate/control_fraction;
ws.surface_rates[p] = - fractions[p] * controls.gas_rate / control_fraction;
}
}
}
@ -1247,8 +1270,8 @@ namespace Opm
}
case Well::ProducerCMode::LRAT:
{
double current_rate = -ws.surface_rates[ pu.phase_pos[Water] ]
- ws.surface_rates[ pu.phase_pos[Oil] ];
Scalar current_rate = - ws.surface_rates[ pu.phase_pos[Water] ]
- ws.surface_rates[ pu.phase_pos[Oil] ];
// or trivial rates or opposite direction we don't just scale the rates
// but use either the potentials or the mobility ratio to initial the well rates
if (current_rate > 0.0) {
@ -1256,8 +1279,8 @@ namespace Opm
ws.surface_rates[p] *= controls.liquid_rate/current_rate;
}
} else {
const std::vector<double> fractions = initialWellRateFractions(simulator, well_state);
double control_fraction = fractions[pu.phase_pos[Water]] + fractions[pu.phase_pos[Oil]];
const std::vector<Scalar> fractions = initialWellRateFractions(simulator, well_state);
const Scalar control_fraction = fractions[pu.phase_pos[Water]] + fractions[pu.phase_pos[Oil]];
if (control_fraction != 0.0) {
for (int p = 0; p<np; ++p) {
ws.surface_rates[p] = - fractions[p] * controls.liquid_rate / control_fraction;
@ -1274,9 +1297,9 @@ namespace Opm
}
case Well::ProducerCMode::RESV:
{
std::vector<double> convert_coeff(this->number_of_phases_, 1.0);
std::vector<Scalar> convert_coeff(this->number_of_phases_, 1.0);
this->rateConverter_.calcCoeff(/*fipreg*/ 0, this->pvtRegionIdx_, ws.surface_rates, convert_coeff);
double total_res_rate = 0.0;
Scalar total_res_rate = 0.0;
for (int p = 0; p<np; ++p) {
total_res_rate -= ws.surface_rates[p] * convert_coeff[p];
}
@ -1288,13 +1311,13 @@ namespace Opm
ws.surface_rates[p] *= controls.resv_rate/total_res_rate;
}
} else {
const std::vector<double> fractions = initialWellRateFractions(simulator, well_state);
const std::vector<Scalar> fractions = initialWellRateFractions(simulator, well_state);
for (int p = 0; p<np; ++p) {
ws.surface_rates[p] = - fractions[p] * controls.resv_rate / convert_coeff[p];
}
}
} else {
std::vector<double> hrates(this->number_of_phases_,0.);
std::vector<Scalar> hrates(this->number_of_phases_,0.);
if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
hrates[pu.phase_pos[Water]] = controls.water_rate;
}
@ -1304,9 +1327,9 @@ namespace Opm
if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
hrates[pu.phase_pos[Gas]] = controls.gas_rate;
}
std::vector<double> hrates_resv(this->number_of_phases_,0.);
std::vector<Scalar> hrates_resv(this->number_of_phases_,0.);
this->rateConverter_.calcReservoirVoidageRates(/*fipreg*/ 0, this->pvtRegionIdx_, hrates, hrates_resv);
double target = std::accumulate(hrates_resv.begin(), hrates_resv.end(), 0.0);
Scalar target = std::accumulate(hrates_resv.begin(), hrates_resv.end(), 0.0);
// or trivial rates or opposite direction we don't just scale the rates
// but use either the potentials or the mobility ratio to initial the well rates
if (total_res_rate > 0.0) {
@ -1314,19 +1337,18 @@ namespace Opm
ws.surface_rates[p] *= target/total_res_rate;
}
} else {
const std::vector<double> fractions = initialWellRateFractions(simulator, well_state);
const std::vector<Scalar> fractions = initialWellRateFractions(simulator, well_state);
for (int p = 0; p<np; ++p) {
ws.surface_rates[p] = - fractions[p] * target / convert_coeff[p];
}
}
}
break;
}
case Well::ProducerCMode::BHP:
{
ws.bhp = controls.bhp_limit;
double total_rate = 0.0;
Scalar total_rate = 0.0;
for (int p = 0; p<np; ++p) {
total_rate -= ws.surface_rates[p];
}
@ -1350,14 +1372,14 @@ namespace Opm
// more sophisticated design might be needed in the future
auto rates = ws.surface_rates;
this->adaptRatesForVFP(rates);
const double bhp = WellBhpThpCalculator(*this).calculateBhpFromThp(
const Scalar bhp = WellBhpThpCalculator(*this).calculateBhpFromThp(
well_state, rates, well, summaryState, this->getRefDensity(), deferred_logger);
ws.bhp = bhp;
ws.thp = this->getTHPConstraint(summaryState);
// if the total rates are negative or zero
// we try to provide a better initial well rate
// using the well potentials
const double total_rate = -std::accumulate(rates.begin(), rates.end(), 0.0);
const Scalar total_rate = -std::accumulate(rates.begin(), rates.end(), 0.0);
if (total_rate <= 0.0) {
for (int p = 0; p < this->number_of_phases_; ++p) {
ws.surface_rates[p] = -ws.well_potentials[p];
@ -1370,14 +1392,14 @@ namespace Opm
{
assert(well.isAvailableForGroupControl());
const auto& group = schedule.getGroup(well.groupName(), this->currentStep());
const double efficiencyFactor = well.getEfficiencyFactor();
double scale = this->getGroupProductionTargetRate(group,
well_state,
group_state,
schedule,
summaryState,
efficiencyFactor,
deferred_logger);
const Scalar efficiencyFactor = well.getEfficiencyFactor();
Scalar scale = this->getGroupProductionTargetRate(group,
well_state,
group_state,
schedule,
summaryState,
efficiencyFactor,
deferred_logger);
// we don't want to scale with zero and get zero rates.
if (scale > 0) {
@ -1394,9 +1416,8 @@ namespace Opm
case Well::ProducerCMode::NONE:
{
OPM_DEFLOG_THROW(std::runtime_error, "Well control must be specified for well " + this->name() , deferred_logger);
}
break;
}
} // end of switch
if (ws.bhp == 0.) {
@ -1406,16 +1427,16 @@ namespace Opm
}
template<typename TypeTag>
std::vector<double>
std::vector<typename WellInterface<TypeTag>::Scalar>
WellInterface<TypeTag>::
initialWellRateFractions(const Simulator& simulator,
const WellState<Scalar>& well_state) const
{
const int np = this->number_of_phases_;
std::vector<double> scaling_factor(np);
std::vector<Scalar> scaling_factor(np);
const auto& ws = well_state.well(this->index_of_well_);
double total_potentials = 0.0;
Scalar total_potentials = 0.0;
for (int p = 0; p<np; ++p) {
total_potentials += ws.well_potentials[p];
}
@ -1427,7 +1448,7 @@ namespace Opm
}
// if we don't have any potentials we weight it using the mobilites
// We only need approximation so we don't bother with the vapporized oil and dissolved gas
double total_tw = 0;
Scalar total_tw = 0;
const int nperf = this->number_of_perforations_;
for (int perf = 0; perf < nperf; ++perf) {
total_tw += this->well_index_[perf];
@ -1436,8 +1457,8 @@ namespace Opm
const int cell_idx = this->well_cells_[perf];
const auto& intQuants = simulator.model().intensiveQuantities(cell_idx, /*timeIdx=*/0);
const auto& fs = intQuants.fluidState();
const double well_tw_fraction = this->well_index_[perf] / total_tw;
double total_mobility = 0.0;
const Scalar well_tw_fraction = this->well_index_[perf] / total_tw;
Scalar total_mobility = 0.0;
for (int p = 0; p < np; ++p) {
int modelPhaseIdx = this->flowPhaseToModelPhaseIdx(p);
total_mobility += fs.invB(modelPhaseIdx).value() * intQuants.mobility(modelPhaseIdx).value();
@ -1463,7 +1484,7 @@ namespace Opm
// if more than one nonzero rate.
auto& ws = well_state.well(this->index_of_well_);
int nonzero_rate_index = -1;
const double floating_point_error_epsilon = 1e-14;
const Scalar floating_point_error_epsilon = 1e-14;
for (int p = 0; p < this->number_of_phases_; ++p) {
if (std::abs(ws.surface_rates[p]) > floating_point_error_epsilon) {
if (nonzero_rate_index == -1) {
@ -1476,7 +1497,7 @@ namespace Opm
}
// Calculate the rates that follow from the current primary variables.
std::vector<double> well_q_s = computeCurrentWellRates(simulator, deferred_logger);
std::vector<Scalar> well_q_s = computeCurrentWellRates(simulator, deferred_logger);
if (nonzero_rate_index == -1) {
// No nonzero rates.
@ -1488,13 +1509,13 @@ namespace Opm
}
// Set the currently-zero phase flows to be nonzero in proportion to well_q_s.
const double initial_nonzero_rate = ws.surface_rates[nonzero_rate_index];
const Scalar initial_nonzero_rate = ws.surface_rates[nonzero_rate_index];
const int comp_idx_nz = this->flowPhaseToModelCompIdx(nonzero_rate_index);
if (std::abs(well_q_s[comp_idx_nz]) > floating_point_error_epsilon) {
for (int p = 0; p < this->number_of_phases_; ++p) {
if (p != nonzero_rate_index) {
const int comp_idx = this->flowPhaseToModelCompIdx(p);
double& rate = ws.surface_rates[p];
Scalar& rate = ws.surface_rates[p];
rate = (initial_nonzero_rate / well_q_s[comp_idx_nz]) * (well_q_s[comp_idx]);
}
}
@ -1502,12 +1523,12 @@ namespace Opm
}
template <typename TypeTag>
std::vector<double>
std::vector<typename WellInterface<TypeTag>::Scalar>
WellInterface<TypeTag>::
wellIndex(const int perf,
const IntensiveQuantities& intQuants,
const double trans_mult,
const SingleWellState<double>& ws) const
const Scalar trans_mult,
const SingleWellState<Scalar>& ws) const
{
// Add a Forchheimer term to the gas phase CTF if the run uses
// either of the WDFAC or the WDFACCOR keywords.
@ -1525,7 +1546,7 @@ namespace Opm
return wi;
}
const double d = this->computeConnectionDFactor(perf, intQuants, ws);
const Scalar d = this->computeConnectionDFactor(perf, intQuants, ws);
if (d < 1.0e-15) {
return wi;
}
@ -1533,43 +1554,43 @@ namespace Opm
// Solve quadratic equations for connection rates satisfying the ipr and the flow-dependent skin.
// If more than one solution, pick the one corresponding to lowest absolute rate (smallest skin).
const auto& connection = this->well_ecl_.getConnections()[ws.perf_data.ecl_index[perf]];
const double Kh = connection.Kh();
const double scaling = 3.141592653589 * Kh * connection.wpimult();
const Scalar Kh = connection.Kh();
const Scalar scaling = 3.141592653589 * Kh * connection.wpimult();
const unsigned gas_comp_idx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
const double connection_pressure = ws.perf_data.pressure[perf];
const double cell_pressure = getValue(intQuants.fluidState().pressure(FluidSystem::gasPhaseIdx));
const double drawdown = cell_pressure - connection_pressure;
const double invB = getValue(intQuants.fluidState().invB(FluidSystem::gasPhaseIdx));
const double mob_g = getValue(intQuants.mobility(FluidSystem::gasPhaseIdx)) * invB;
const double a = d;
const double b = 2*scaling/wi[gas_comp_idx];
const double c = -2*scaling*mob_g*drawdown;
const Scalar connection_pressure = ws.perf_data.pressure[perf];
const Scalar cell_pressure = getValue(intQuants.fluidState().pressure(FluidSystem::gasPhaseIdx));
const Scalar drawdown = cell_pressure - connection_pressure;
const Scalar invB = getValue(intQuants.fluidState().invB(FluidSystem::gasPhaseIdx));
const Scalar mob_g = getValue(intQuants.mobility(FluidSystem::gasPhaseIdx)) * invB;
const Scalar a = d;
const Scalar b = 2*scaling/wi[gas_comp_idx];
const Scalar c = -2*scaling*mob_g*drawdown;
double consistent_Q = -1.0e20;
Scalar consistent_Q = -1.0e20;
// Find and check negative solutions (a --> -a)
const double r2n = b*b + 4*a*c;
const Scalar r2n = b*b + 4*a*c;
if (r2n >= 0) {
const double rn = std::sqrt(r2n);
const double xn1 = (b-rn)*0.5/a;
const Scalar rn = std::sqrt(r2n);
const Scalar xn1 = (b-rn)*0.5/a;
if (xn1 <= 0) {
consistent_Q = xn1;
}
const double xn2 = (b+rn)*0.5/a;
const Scalar xn2 = (b+rn)*0.5/a;
if (xn2 <= 0 && xn2 > consistent_Q) {
consistent_Q = xn2;
}
}
// Find and check positive solutions
consistent_Q *= -1;
const double r2p = b*b - 4*a*c;
const Scalar r2p = b*b - 4*a*c;
if (r2p >= 0) {
const double rp = std::sqrt(r2p);
const double xp1 = (rp-b)*0.5/a;
const Scalar rp = std::sqrt(r2p);
const Scalar xp1 = (rp-b)*0.5/a;
if (xp1 > 0 && xp1 < consistent_Q) {
consistent_Q = xp1;
}
const double xp2 = -(rp+b)*0.5/a;
const Scalar xp2 = -(rp+b)*0.5/a;
if (xp2 > 0 && xp2 < consistent_Q) {
consistent_Q = xp2;
}
@ -1583,7 +1604,7 @@ namespace Opm
void
WellInterface<TypeTag>::
updateConnectionDFactor(const Simulator& simulator,
SingleWellState<double>& ws) const
SingleWellState<Scalar>& ws) const
{
if (! this->well_ecl_.getWDFAC().useDFactor()) {
return;
@ -1600,11 +1621,11 @@ namespace Opm
}
template <typename TypeTag>
double
typename WellInterface<TypeTag>::Scalar
WellInterface<TypeTag>::
computeConnectionDFactor(const int perf,
const IntensiveQuantities& intQuants,
const SingleWellState<double>& ws) const
const SingleWellState<Scalar>& ws) const
{
auto rhoGS = [regIdx = this->pvtRegionIdx()]() {
return FluidSystem::referenceDensity(FluidSystem::gasPhaseIdx, regIdx);
@ -1622,7 +1643,7 @@ namespace Opm
// Note that rv here is from grid block with typically
// p_block > connection_pressure
// so we may very well have rv > rv_sat
const double rv_sat = gasPvt.saturatedOilVaporizationFactor
const Scalar rv_sat = gasPvt.saturatedOilVaporizationFactor
(regIdx, temperature, connection_pressure);
if (! (rv < rv_sat)) {
@ -1645,7 +1666,7 @@ namespace Opm
void
WellInterface<TypeTag>::
updateConnectionTransmissibilityFactor(const Simulator& simulator,
SingleWellState<double>& ws) const
SingleWellState<Scalar>& ws) const
{
auto connCF = [&connIx = std::as_const(ws.perf_data.ecl_index),
&conns = this->well_ecl_.getConnections()]
@ -1772,7 +1793,7 @@ namespace Opm
auto bhp_at_thp_limit = computeBhpAtThpLimitProdWithAlq(
simulator, summary_state, this->getALQ(well_state), deferred_logger);
if (bhp_at_thp_limit) {
std::vector<double> rates(this->number_of_phases_, 0.0);
std::vector<Scalar> rates(this->number_of_phases_, 0.0);
if (thp_update_iterations) {
computeWellRatesWithBhpIterations(simulator, *bhp_at_thp_limit,
rates, deferred_logger);
@ -1794,9 +1815,9 @@ namespace Opm
void
WellInterface<TypeTag>::
computeConnLevelProdInd(const FluidState& fs,
const std::function<double(const double)>& connPICalc,
const std::function<Scalar(const Scalar)>& connPICalc,
const std::vector<Scalar>& mobility,
double* connPI) const
Scalar* connPI) const
{
const auto& pu = this->phaseUsage();
const int np = this->number_of_phases_;
@ -1830,9 +1851,9 @@ namespace Opm
WellInterface<TypeTag>::
computeConnLevelInjInd(const FluidState& fs,
const Phase preferred_phase,
const std::function<double(const double)>& connIICalc,
const std::function<Scalar(const Scalar)>& connIICalc,
const std::vector<Scalar>& mobility,
double* connII,
Scalar* connII,
DeferredLogger& deferred_logger) const
{
// Assumes single phase injection

6
opm/simulators/wells/WellTest.hpp
View File

@ -34,14 +34,14 @@ class DeferredLogger;
struct PhaseUsage;
template<class Scalar> class SingleWellState;
class WellEconProductionLimits;
class WellInterfaceGeneric;
template<class Scalar> class WellInterfaceGeneric;
class WellTestState;
//! \brief Class for performing well tests.
class WellTest {
public:
//! \brief Constructor sets reference to well.
WellTest(const WellInterfaceGeneric& well) : well_(well) {}
WellTest(const WellInterfaceGeneric<double>& well) : well_(well) {}
void updateWellTestStateEconomic(const SingleWellState<double>& ws,
const double simulation_time,
@ -95,7 +95,7 @@ private:
DeferredLogger& deferred_logger) const;
const WellInterfaceGeneric& well_; //!< Reference to well interface
const WellInterfaceGeneric<double>& well_; //!< Reference to well interface
};
}