diff --git a/opm/simulators/wells/BlackoilWellModelGeneric.cpp b/opm/simulators/wells/BlackoilWellModelGeneric.cpp index fb4ef6e07..1e6a25fe2 100644 --- a/opm/simulators/wells/BlackoilWellModelGeneric.cpp +++ b/opm/simulators/wells/BlackoilWellModelGeneric.cpp @@ -1334,14 +1334,14 @@ calculateEfficiencyFactors(const int reportStepIdx) } template -WellInterfaceGeneric* +WellInterfaceGeneric* BlackoilWellModelGeneric:: 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* elem)->bool { return elem->name() == well_name; }); diff --git a/opm/simulators/wells/BlackoilWellModelGeneric.hpp b/opm/simulators/wells/BlackoilWellModelGeneric.hpp index d3ad45c85..45875a5a7 100644 --- a/opm/simulators/wells/BlackoilWellModelGeneric.hpp +++ b/opm/simulators/wells/BlackoilWellModelGeneric.hpp @@ -64,7 +64,7 @@ namespace Opm { struct SimulatorUpdate; class SummaryConfig; class VFPProperties; - class WellInterfaceGeneric; + template class WellInterfaceGeneric; template class WellState; } // namespace Opm @@ -84,7 +84,7 @@ class BlackoilWellModelGeneric public: // --------- Types --------- using GLiftOptWells = std::map>>; - using GLiftProdWells = std::map; + using GLiftProdWells = std::map*>; using GLiftWellStateMap = std::map>>; BlackoilWellModelGeneric(Schedule& schedule, @@ -115,7 +115,7 @@ public: const Schedule& schedule() const { return schedule_; } const PhaseUsage& phaseUsage() const { return phase_usage_; } const GroupState& groupState() const { return this->active_wgstate_.group_state; } - std::vector genericWells() const + std::vector*> genericWells() const { return {well_container_generic_.begin(), well_container_generic_.end()}; } /* @@ -542,7 +542,7 @@ protected: std::function not_on_process_{}; // a vector of all the wells. - std::vector well_container_generic_{}; + std::vector*> well_container_generic_{}; std::vector local_shut_wells_{}; @@ -589,7 +589,7 @@ protected: std::map> closed_offending_wells_; private: - WellInterfaceGeneric* getGenWell(const std::string& well_name); + WellInterfaceGeneric* getGenWell(const std::string& well_name); }; diff --git a/opm/simulators/wells/BlackoilWellModelGuideRates.cpp b/opm/simulators/wells/BlackoilWellModelGuideRates.cpp index 6188181a0..3b9cccb52 100644 --- a/opm/simulators/wells/BlackoilWellModelGuideRates.cpp +++ b/opm/simulators/wells/BlackoilWellModelGuideRates.cpp @@ -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* 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* well) { return events.hasEvent(well->name(), effective_events_mask); }); diff --git a/opm/simulators/wells/GasLiftSingleWell.hpp b/opm/simulators/wells/GasLiftSingleWell.hpp index 5fab8c89b..4ce75dcef 100644 --- a/opm/simulators/wells/GasLiftSingleWell.hpp +++ b/opm/simulators/wells/GasLiftSingleWell.hpp @@ -51,7 +51,7 @@ public: const Parallel::Communication& comm, bool glift_debug); - const WellInterfaceGeneric& getWell() const override { return well_; } + const WellInterfaceGeneric& getWell() const override { return well_; } private: std::optional diff --git a/opm/simulators/wells/GasLiftSingleWellGeneric.hpp b/opm/simulators/wells/GasLiftSingleWellGeneric.hpp index 6fac11ad2..8341b6c18 100644 --- a/opm/simulators/wells/GasLiftSingleWellGeneric.hpp +++ b/opm/simulators/wells/GasLiftSingleWellGeneric.hpp @@ -40,7 +40,7 @@ class GasLiftWell; template class GasLiftWellState; class Schedule; class SummaryState; -class WellInterfaceGeneric; +template class WellInterfaceGeneric; template class WellState; template class GroupState; @@ -107,7 +107,7 @@ public: std::unique_ptr> runOptimize(const int iteration_idx); - virtual const WellInterfaceGeneric& getWell() const = 0; + virtual const WellInterfaceGeneric& getWell() const = 0; protected: GasLiftSingleWellGeneric(DeferredLogger& deferred_logger, diff --git a/opm/simulators/wells/GasLiftSingleWell_impl.hpp b/opm/simulators/wells/GasLiftSingleWell_impl.hpp index adc83bef1..bdab93459 100644 --- a/opm/simulators/wells/GasLiftSingleWell_impl.hpp +++ b/opm/simulators/wells/GasLiftSingleWell_impl.hpp @@ -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) { diff --git a/opm/simulators/wells/GasLiftStage2.cpp b/opm/simulators/wells/GasLiftStage2.cpp index 36d5331ef..74b7faff4 100644 --- a/opm/simulators/wells/GasLiftStage2.cpp +++ b/opm/simulators/wells/GasLiftStage2.cpp @@ -1118,7 +1118,7 @@ computeDelta(const std::string& well_name) const GradInfo& gi = this->parent.dec_grads_.at(well_name); GasLiftWellState& 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& 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(); diff --git a/opm/simulators/wells/GasLiftStage2.hpp b/opm/simulators/wells/GasLiftStage2.hpp index 0611147cb..f6b86df90 100644 --- a/opm/simulators/wells/GasLiftStage2.hpp +++ b/opm/simulators/wells/GasLiftStage2.hpp @@ -38,7 +38,7 @@ template class GasLiftWellState; class Group; template class GroupState; class Schedule; -class WellInterfaceGeneric; +template class WellInterfaceGeneric; template class WellState; template @@ -46,7 +46,7 @@ class GasLiftStage2 : public GasLiftCommon { using GasLiftSingleWell = GasLiftSingleWellGeneric; using GLiftOptWells = std::map>; - using GLiftProdWells = std::map; + using GLiftProdWells = std::map*>; using GLiftWellStateMap = std::map>>; using GradPair = std::pair; using GradPairItr = typename std::vector::iterator; diff --git a/opm/simulators/wells/MultisegmentWellEquations.cpp b/opm/simulators/wells/MultisegmentWellEquations.cpp index dd0644bd6..0f86ff51b 100644 --- a/opm/simulators/wells/MultisegmentWellEquations.cpp +++ b/opm/simulators/wells/MultisegmentWellEquations.cpp @@ -311,7 +311,7 @@ extractCPRPressureMatrix(PressureMatrix& jacobian, const BVector& weights, const int pressureVarIndex, const bool /*use_well_weights*/, - const WellInterfaceGeneric& well, + const WellInterfaceGeneric& well, const int seg_pressure_var_ind, const WellState& well_state) const { @@ -395,7 +395,7 @@ template void MultisegmentWellEquations:: \ const MultisegmentWellEquations::BVector&, \ const int, \ const bool, \ - const WellInterfaceGeneric&, \ + const WellInterfaceGeneric&, \ const int, \ const WellState&) const; diff --git a/opm/simulators/wells/MultisegmentWellEquations.hpp b/opm/simulators/wells/MultisegmentWellEquations.hpp index 8a06655a3..58289e48e 100644 --- a/opm/simulators/wells/MultisegmentWellEquations.hpp +++ b/opm/simulators/wells/MultisegmentWellEquations.hpp @@ -41,7 +41,7 @@ template class MultisegmentWellGeneric; #if COMPILE_BDA_BRIDGE class WellContributions; #endif -class WellInterfaceGeneric; +template class WellInterfaceGeneric; template class WellState; template @@ -118,7 +118,7 @@ public: const BVector& weights, const int pressureVarIndex, const bool /*use_well_weights*/, - const WellInterfaceGeneric& well, + const WellInterfaceGeneric& well, const int seg_pressure_var_ind, const WellState& well_state) const; diff --git a/opm/simulators/wells/MultisegmentWellGeneric.cpp b/opm/simulators/wells/MultisegmentWellGeneric.cpp index 338070101..994bd0d0b 100644 --- a/opm/simulators/wells/MultisegmentWellGeneric.cpp +++ b/opm/simulators/wells/MultisegmentWellGeneric.cpp @@ -34,18 +34,15 @@ #include #include -#include #include -#include #include -namespace Opm -{ +namespace Opm { template MultisegmentWellGeneric:: -MultisegmentWellGeneric(WellInterfaceGeneric& baseif) +MultisegmentWellGeneric(WellInterfaceGeneric& baseif) : baseif_(baseif) { } diff --git a/opm/simulators/wells/MultisegmentWellGeneric.hpp b/opm/simulators/wells/MultisegmentWellGeneric.hpp index 678067e6c..a2509feb5 100644 --- a/opm/simulators/wells/MultisegmentWellGeneric.hpp +++ b/opm/simulators/wells/MultisegmentWellGeneric.hpp @@ -32,7 +32,7 @@ namespace Opm class DeferredLogger; class SummaryState; -class WellInterfaceGeneric; +template class WellInterfaceGeneric; enum class WellSegmentCompPressureDrop; class WellSegments; template class WellState; @@ -52,7 +52,7 @@ public: int numberOfSegments() const; protected: - MultisegmentWellGeneric(WellInterfaceGeneric& baseif); + MultisegmentWellGeneric(WellInterfaceGeneric& baseif); // scale the segment rates and pressure based on well rates and bhp void scaleSegmentRatesWithWellRates(const std::vector>& segment_inlets, @@ -75,7 +75,7 @@ protected: const double density, const std::vector& seg_dp) const; - const WellInterfaceGeneric& baseif_; + const WellInterfaceGeneric& baseif_; }; } diff --git a/opm/simulators/wells/MultisegmentWellSegments.cpp b/opm/simulators/wells/MultisegmentWellSegments.cpp index 95d2a893f..dcf208911 100644 --- a/opm/simulators/wells/MultisegmentWellSegments.cpp +++ b/opm/simulators/wells/MultisegmentWellSegments.cpp @@ -45,7 +45,6 @@ #include -#include #include #include #include @@ -61,7 +60,7 @@ namespace Opm template MultisegmentWellSegments:: MultisegmentWellSegments(const int numSegments, - WellInterfaceGeneric& well) + WellInterfaceGeneric& well) : perforations_(numSegments) , perforation_depth_diffs_(well.numPerfs(), 0.0) , inlets_(well.wellEcl().getSegments().size()) diff --git a/opm/simulators/wells/MultisegmentWellSegments.hpp b/opm/simulators/wells/MultisegmentWellSegments.hpp index 4ae2d75c9..eaa51b397 100644 --- a/opm/simulators/wells/MultisegmentWellSegments.hpp +++ b/opm/simulators/wells/MultisegmentWellSegments.hpp @@ -33,7 +33,7 @@ namespace Opm { struct PhaseUsage; template class SegmentState; class UnitSystem; - class WellInterfaceGeneric; + template class WellInterfaceGeneric; class SummaryState; } // namespace Opm @@ -49,7 +49,7 @@ class MultisegmentWellSegments public: MultisegmentWellSegments(const int numSegments, - WellInterfaceGeneric& well); + WellInterfaceGeneric& well); void computeFluidProperties(const EvalWell& temperature, const EvalWell& saltConcentration, @@ -172,7 +172,7 @@ private: std::vector> phase_fractions_; std::vector> phase_viscosities_; - WellInterfaceGeneric& well_; + WellInterfaceGeneric& well_; void copyPhaseDensities(const unsigned phaseIdx, const std::size_t stride, diff --git a/opm/simulators/wells/MultisegmentWell_impl.hpp b/opm/simulators/wells/MultisegmentWell_impl.hpp index a2e936bb7..d62f35e03 100644 --- a/opm/simulators/wells/MultisegmentWell_impl.hpp +++ b/opm/simulators/wells/MultisegmentWell_impl.hpp @@ -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::computeWellPotentials(well_potentials, well_state); if (!compute_potential) { return; diff --git a/opm/simulators/wells/StandardWellEquations.cpp b/opm/simulators/wells/StandardWellEquations.cpp index f819d4c68..86742d3ca 100644 --- a/opm/simulators/wells/StandardWellEquations.cpp +++ b/opm/simulators/wells/StandardWellEquations.cpp @@ -291,7 +291,7 @@ extractCPRPressureMatrix(PressureMatrix& jacobian, const BVector& weights, const int pressureVarIndex, const bool use_well_weights, - const WellInterfaceGeneric& well, + const WellInterfaceGeneric& well, const int bhp_var_index, const WellState& well_state) const { @@ -413,7 +413,7 @@ template void StandardWellEquations:: \ const typename StandardWellEquations::BVector&, \ const int, \ const bool, \ - const WellInterfaceGeneric&, \ + const WellInterfaceGeneric&, \ const int, \ const WellState&) const; diff --git a/opm/simulators/wells/StandardWellEquations.hpp b/opm/simulators/wells/StandardWellEquations.hpp index 18935eae9..8636c376d 100644 --- a/opm/simulators/wells/StandardWellEquations.hpp +++ b/opm/simulators/wells/StandardWellEquations.hpp @@ -39,7 +39,7 @@ template class StandardWellEquationAccess; #if COMPILE_BDA_BRIDGE class WellContributions; #endif -class WellInterfaceGeneric; +template class WellInterfaceGeneric; template class WellState; template @@ -115,7 +115,7 @@ public: const BVector& weights, const int pressureVarIndex, const bool use_well_weights, - const WellInterfaceGeneric& well, + const WellInterfaceGeneric& well, const int bhp_var_index, const WellState& well_state) const; diff --git a/opm/simulators/wells/StandardWell_impl.hpp b/opm/simulators/wells/StandardWell_impl.hpp index ef0c7aea9..5331db438 100644 --- a/opm/simulators/wells/StandardWell_impl.hpp +++ b/opm/simulators/wells/StandardWell_impl.hpp @@ -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::computeWellPotentials(well_potentials, well_state); if (!compute_potential) { return; diff --git a/opm/simulators/wells/WellBhpThpCalculator.hpp b/opm/simulators/wells/WellBhpThpCalculator.hpp index 2b8fb09cb..fdf0502ca 100644 --- a/opm/simulators/wells/WellBhpThpCalculator.hpp +++ b/opm/simulators/wells/WellBhpThpCalculator.hpp @@ -35,14 +35,14 @@ namespace Opm class DeferredLogger; class SummaryState; class Well; -class WellInterfaceGeneric; +template class WellInterfaceGeneric; template 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& 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& well_; //!< Reference to well interface }; } diff --git a/opm/simulators/wells/WellConstraints.hpp b/opm/simulators/wells/WellConstraints.hpp index b33ac748b..5bec0b6b9 100644 --- a/opm/simulators/wells/WellConstraints.hpp +++ b/opm/simulators/wells/WellConstraints.hpp @@ -39,7 +39,7 @@ using RegionId = int; class Rates; template class SingleWellState; class SummaryState; -class WellInterfaceGeneric; +template 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& well) : well_(well) {} using RateConvFunc = std::function&, @@ -78,7 +78,7 @@ private: DeferredLogger& deferred_logger, const std::optional& prod_controls = std::nullopt) const; - const WellInterfaceGeneric& well_; //!< Reference to well interface + const WellInterfaceGeneric& well_; //!< Reference to well interface }; } diff --git a/opm/simulators/wells/WellConvergence.hpp b/opm/simulators/wells/WellConvergence.hpp index 27beadac4..e510f2631 100644 --- a/opm/simulators/wells/WellConvergence.hpp +++ b/opm/simulators/wells/WellConvergence.hpp @@ -30,13 +30,13 @@ namespace Opm class ConvergenceReport; class DeferredLogger; -class WellInterfaceGeneric; +template class WellInterfaceGeneric; template class WellState; class WellConvergence { public: - WellConvergence(const WellInterfaceGeneric& well) + WellConvergence(const WellInterfaceGeneric& well) : well_(well) {} @@ -62,7 +62,7 @@ public: ConvergenceReport& report) const; private: - const WellInterfaceGeneric& well_; + const WellInterfaceGeneric& well_; }; } diff --git a/opm/simulators/wells/WellFilterCake.cpp b/opm/simulators/wells/WellFilterCake.cpp index c1d8dff9a..9ae8f4e2c 100644 --- a/opm/simulators/wells/WellFilterCake.cpp +++ b/opm/simulators/wells/WellFilterCake.cpp @@ -36,7 +36,7 @@ namespace Opm { template void WellFilterCake:: -updateFiltrationParticleVolume(const WellInterfaceGeneric& well, +updateFiltrationParticleVolume(const WellInterfaceGeneric& well, const double dt, const Scalar conc, const std::size_t water_index, @@ -78,7 +78,7 @@ updateFiltrationParticleVolume(const WellInterfaceGeneric& well, template void WellFilterCake:: -updateInjFCMult(const WellInterfaceGeneric& well, +updateInjFCMult(const WellInterfaceGeneric& well, WellState& well_state, DeferredLogger& deferred_logger) { diff --git a/opm/simulators/wells/WellFilterCake.hpp b/opm/simulators/wells/WellFilterCake.hpp index 863fc7b7f..b9e8d37ae 100644 --- a/opm/simulators/wells/WellFilterCake.hpp +++ b/opm/simulators/wells/WellFilterCake.hpp @@ -26,7 +26,7 @@ namespace Opm { class DeferredLogger; -class WellInterfaceGeneric; +template class WellInterfaceGeneric; template 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& well, const double dt, const Scalar conc, const std::size_t water_index, WellState& well_state); //! \brief Update the multiplier for well transmissbility due to cake filtration. - void updateInjFCMult(const WellInterfaceGeneric& well, + void updateInjFCMult(const WellInterfaceGeneric& well, WellState& well_state, DeferredLogger& deferred_logger); diff --git a/opm/simulators/wells/WellGroupConstraints.hpp b/opm/simulators/wells/WellGroupConstraints.hpp index 0176cded9..5e9a3d272 100644 --- a/opm/simulators/wells/WellGroupConstraints.hpp +++ b/opm/simulators/wells/WellGroupConstraints.hpp @@ -40,16 +40,19 @@ enum class InjectorType; using RegionId = int; class Schedule; class SummaryState; -class WellInterfaceGeneric; +template class WellInterfaceGeneric; template 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& well) : well_(well) {} - using RateConvFunc = std::function&, std::vector&)>; + using RateConvFunc = std::function&, + std::vector&)>; bool checkGroupConstraints(WellState& well_state, const GroupState& group_state, @@ -79,7 +82,7 @@ private: const RateConvFunc& rateConverter, DeferredLogger& deferred_logger) const; - const WellInterfaceGeneric& well_; //!< Reference to well interface + const WellInterfaceGeneric& well_; //!< Reference to well interface }; } diff --git a/opm/simulators/wells/WellGroupControls.hpp b/opm/simulators/wells/WellGroupControls.hpp index 9563eda37..899854ee3 100644 --- a/opm/simulators/wells/WellGroupControls.hpp +++ b/opm/simulators/wells/WellGroupControls.hpp @@ -39,14 +39,14 @@ enum class InjectorType; using RegionId = int; class Schedule; class SummaryState; -class WellInterfaceGeneric; +template class WellInterfaceGeneric; template 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& well) : well_(well) {} using RateConvFunc = std::function&, std::vector&)>; @@ -98,7 +98,7 @@ public: DeferredLogger& deferred_logger) const; private: - const WellInterfaceGeneric& well_; //!< Reference to well interface + const WellInterfaceGeneric& well_; //!< Reference to well interface }; } diff --git a/opm/simulators/wells/WellInterface.hpp b/opm/simulators/wells/WellInterface.hpp index eaa22bda9..25018ed43 100644 --- a/opm/simulators/wells/WellInterface.hpp +++ b/opm/simulators/wells/WellInterface.hpp @@ -98,7 +98,7 @@ public: using MatrixBlockType = Dune::FieldMatrix; using Eval = typename Base::Eval; using BVector = Dune::BlockVector; - using PressureMatrix = Dune::BCRSMatrix>; + using PressureMatrix = Dune::BCRSMatrix>; using RateConverterType = typename WellInterfaceFluidSystem::RateConverterType; @@ -148,17 +148,17 @@ public: virtual ~WellInterface() = default; virtual void init(const PhaseUsage* phase_usage_arg, - const std::vector& depth_arg, - const double gravity_arg, + const std::vector& depth_arg, + const Scalar gravity_arg, const int num_cells, - const std::vector< Scalar >& B_avg, + const std::vector& B_avg, const bool changed_to_open_this_step); virtual void initPrimaryVariablesEvaluation() = 0; virtual ConvergenceReport getWellConvergence(const SummaryState& summary_state, const WellState& well_state, - const std::vector& B_avg, + const std::vector& 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& well_flux, - DeferredLogger& deferred_logger - ) const = 0; + virtual void computeWellRatesWithBhp(const Simulator& ebosSimulator, + const Scalar& bhp, + std::vector& well_flux, + DeferredLogger& deferred_logger) const = 0; - virtual std::optional computeBhpAtThpLimitProdWithAlq( - const Simulator& simulator, - const SummaryState& summary_state, - const double alq_value, - DeferredLogger& deferred_logger - ) const = 0; + virtual std::optional + 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& well_state, - std::vector& well_potentials, + std::vector& 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& well_flux, + std::vector& well_flux, DeferredLogger& deferred_logger) const = 0; bool updateWellStateWithTHPTargetProd(const Simulator& simulator, @@ -245,7 +242,7 @@ public: const GroupState& 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& 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& well_state, DeferredLogger& deferred_logger); - bool gliftBeginTimeStepWellTestIterateWellEquations( - const Simulator& simulator, - const double dt, - WellState& well_state, - const GroupState& group_state, - DeferredLogger& deferred_logger); + bool gliftBeginTimeStepWellTestIterateWellEquations(const Simulator& ebos_simulator, + const double dt, + WellState& well_state, + const GroupState& group_state, + DeferredLogger& deferred_logger); void gliftBeginTimeStepWellTestUpdateALQ(const Simulator& simulator, WellState& well_state, @@ -325,8 +321,9 @@ public: /// Compute well rates based on current reservoir conditions and well variables. /// Used in updateWellStateRates(). - virtual std::vector computeCurrentWellRates(const Simulator& simulator, - DeferredLogger& deferred_logger) const = 0; + virtual std::vector + 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 getPrimaryVars() const + virtual std::vector getPrimaryVars() const { return {}; } - virtual int setPrimaryVars(std::vector::const_iterator) + virtual int setPrimaryVars(typename std::vector::const_iterator) { return 0; } - std::vector wellIndex(const int perf, + std::vector wellIndex(const int perf, const IntensiveQuantities& intQuants, - const double trans_mult, - const SingleWellState& ws) const; + const Scalar trans_mult, + const SingleWellState& ws) const; void updateConnectionDFactor(const Simulator& simulator, - SingleWellState& ws) const; + SingleWellState& ws) const; void updateConnectionTransmissibilityFactor(const Simulator& simulator, - SingleWellState& ws) const; - + SingleWellState& ws) const; protected: // simulation parameters const ModelParameters& param_; std::vector connectionRates_; - std::vector< Scalar > B_avg_; + std::vector 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& compFrac() const; + const std::vector& compFrac() const; - std::vector initialWellRateFractions(const Simulator& simulator, - const WellState& well_state) const; + std::vector + initialWellRateFractions(const Simulator& ebosSimulator, + const WellState& well_state) const; // check whether the well is operable under BHP limit with current reservoir condition virtual void checkOperabilityUnderBHPLimit(const WellState& well_state, @@ -453,15 +445,16 @@ protected: const GroupState& group_state, DeferredLogger& deferred_logger); - std::optional estimateOperableBhp(const Simulator& simulator, - const double dt, - WellState& well_state, - const SummaryState& summary_state, - DeferredLogger& deferred_logger); + std::optional + estimateOperableBhp(const Simulator& ebos_simulator, + const double dt, + WellState& well_state, + const SummaryState& summary_state, + DeferredLogger& deferred_logger); bool solveWellWithBhp(const Simulator& simulator, const double dt, - const double bhp, + const Scalar bhp, WellState& well_state, DeferredLogger& deferred_logger); @@ -486,21 +479,20 @@ protected: [[maybe_unused]] DeferredLogger& deferred_logger) const; void computeConnLevelProdInd(const FluidState& fs, - const std::function& connPICalc, + const std::function& connPICalc, const std::vector& mobility, - double* connPI) const; + Scalar* connPI) const; void computeConnLevelInjInd(const FluidState& fs, const Phase preferred_phase, - const std::function& connIICalc, + const std::function& connIICalc, const std::vector& mobility, - double* connII, + Scalar* connII, DeferredLogger& deferred_logger) const; - double computeConnectionDFactor(const int perf, + Scalar computeConnectionDFactor(const int perf, const IntensiveQuantities& intQuants, - const SingleWellState& ws) const; - + const SingleWellState& ws) const; }; } // namespace Opm diff --git a/opm/simulators/wells/WellInterfaceFluidSystem.cpp b/opm/simulators/wells/WellInterfaceFluidSystem.cpp index 6e90e9256..51aeb45db 100644 --- a/opm/simulators/wells/WellInterfaceFluidSystem.cpp +++ b/opm/simulators/wells/WellInterfaceFluidSystem.cpp @@ -54,9 +54,9 @@ WellInterfaceFluidSystem(const Well& well, const int num_phases, const int index_of_well, const std::vector& perf_data) - : WellInterfaceGeneric(well, parallel_well_info, time_step, - pvtRegionIdx, num_components, num_phases, - index_of_well, perf_data) + : WellInterfaceGeneric(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 void WellInterfaceFluidSystem:: -calculateReservoirRates(SingleWellState& ws) const +calculateReservoirRates(SingleWellState& 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& 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 + const auto surface_rates_perf = std::vector { surf_perf_rates.begin() + (i + 0)*np , surf_perf_rates.begin() + (i + 1)*np }; - std::vector voidage_rates_perf(np, 0.0); + std::vector voidage_rates_perf(np, 0.0); this->rateConverter_ .calcReservoirVoidageRates(fipreg, this->pvtRegionIdx_, @@ -101,11 +101,11 @@ calculateReservoirRates(SingleWellState& 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& 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 + const auto surface_rates_perf = std::vector { 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 voidage_rates_perf(np, 0.0); + std::vector voidage_rates_perf(np, 0.0); this->rateConverter_ .calcReservoirVoidageRates(this->pvtRegionIdx_, pressure, @@ -153,7 +153,7 @@ calculateReservoirRates(SingleWellState& ws) const template bool WellInterfaceFluidSystem:: -checkIndividualConstraints(SingleWellState& ws, +checkIndividualConstraints(SingleWellState& ws, const SummaryState& summaryState, DeferredLogger& deferred_logger, const std::optional& inj_controls, @@ -161,8 +161,8 @@ checkIndividualConstraints(SingleWellState& ws, { auto rRates = [this](const int fipreg, const int pvtRegion, - const std::vector& surface_rates, - std::vector& voidage_rates) + const std::vector& surface_rates, + std::vector& voidage_rates) { return rateConverter_.calcReservoirVoidageRates(fipreg, pvtRegion, surface_rates, voidage_rates); @@ -177,23 +177,25 @@ checkIndividualConstraints(SingleWellState& ws, template bool WellInterfaceFluidSystem:: -checkGroupConstraints(WellState& well_state, - const GroupState& group_state, +checkGroupConstraints(WellState& well_state, + const GroupState& 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& prod_gname, std::vector& coeff) + auto rCoeff = [this, &group_state](const RegionId id, + const int region, + const std::optional& prod_gname, + std::vector& 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& well_state, template bool WellInterfaceFluidSystem:: -checkConstraints(WellState& well_state, - const GroupState& group_state, +checkConstraints(WellState& well_state, + const GroupState& 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 -std::optional +std::optional::Scalar> WellInterfaceFluidSystem:: getGroupInjectionTargetRate(const Group& group, - const WellState& well_state, - const GroupState& group_state, + const WellState& well_state, + const GroupState& 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& prod_gname, std::vector& coeff) + auto rCoeff = [this, &group_state](const RegionId id, const int region, + const std::optional& prod_gname, + std::vector& 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 -double +typename WellInterfaceFluidSystem::Scalar WellInterfaceFluidSystem:: getGroupProductionTargetRate(const Group& group, - const WellState& well_state, - const GroupState& group_state, + const WellState& well_state, + const GroupState& 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& prod_gname, std::vector& coeff) + auto rCoeff = [this, &group_state](const RegionId id, const int region, + const std::optional& prod_gname, + std::vector& 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); } diff --git a/opm/simulators/wells/WellInterfaceFluidSystem.hpp b/opm/simulators/wells/WellInterfaceFluidSystem.hpp index dea30feff..4949066dd 100644 --- a/opm/simulators/wells/WellInterfaceFluidSystem.hpp +++ b/opm/simulators/wells/WellInterfaceFluidSystem.hpp @@ -44,7 +44,8 @@ template class SingleWellState; template class WellState; template -class WellInterfaceFluidSystem : public WellInterfaceGeneric { +class WellInterfaceFluidSystem : public WellInterfaceGeneric +{ protected: using RateConverterType = RateConverter:: SurfaceToReservoirVoidage>; @@ -52,6 +53,8 @@ protected: static constexpr int INVALIDCOMPLETION = std::numeric_limits::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& perf_data); // updating the voidage rates in well_state when requested - void calculateReservoirRates(SingleWellState& ws) const; + void calculateReservoirRates(SingleWellState& ws) const; - bool checkIndividualConstraints(SingleWellState& ws, + bool checkIndividualConstraints(SingleWellState& ws, const SummaryState& summaryState, DeferredLogger& deferred_logger, const std::optional& inj_controls = std::nullopt, const std::optional& prod_controls = std::nullopt) const; - bool checkGroupConstraints(WellState& well_state, - const GroupState& group_state, + bool checkGroupConstraints(WellState& well_state, + const GroupState& group_state, const Schedule& schedule, const SummaryState& summaryState, DeferredLogger& deferred_logger) const; - bool checkConstraints(WellState& well_state, - const GroupState& group_state, + bool checkConstraints(WellState& well_state, + const GroupState& group_state, const Schedule& schedule, const SummaryState& summaryState, DeferredLogger& deferred_logger) const; - std::optional + std::optional getGroupInjectionTargetRate(const Group& group, - const WellState& well_state, - const GroupState& group_state, + const WellState& well_state, + const GroupState& 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& well_state, - const GroupState& group_state, + const WellState& well_state, + const GroupState& 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 diff --git a/opm/simulators/wells/WellInterfaceGeneric.cpp b/opm/simulators/wells/WellInterfaceGeneric.cpp index 5fbc34609..2cf2faec3 100644 --- a/opm/simulators/wells/WellInterfaceGeneric.cpp +++ b/opm/simulators/wells/WellInterfaceGeneric.cpp @@ -49,17 +49,18 @@ #include #include -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& perf_data) +template +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& 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& rates) const +template +void WellInterfaceGeneric:: +adaptRatesForVFP(std::vector& rates) const { const auto& pu = this->phaseUsage(); if (pu.num_phases == 2) { @@ -139,73 +142,91 @@ void WellInterfaceGeneric::adaptRatesForVFP(std::vector& rates) const } } -const std::vector& WellInterfaceGeneric::perforationData() const +template +const std::vector& +WellInterfaceGeneric::perforationData() const { return *perf_data_; } -const std::string& WellInterfaceGeneric::name() const +template +const std::string& +WellInterfaceGeneric::name() const { return well_ecl_.name(); } -bool WellInterfaceGeneric::isInjector() const +template +bool WellInterfaceGeneric::isInjector() const { return well_ecl_.isInjector(); } -bool WellInterfaceGeneric::isProducer() const +template +bool WellInterfaceGeneric::isProducer() const { return well_ecl_.isProducer(); } -int WellInterfaceGeneric::indexOfWell() const +template +int WellInterfaceGeneric::indexOfWell() const { return index_of_well_; } -bool WellInterfaceGeneric::getAllowCrossFlow() const +template +bool WellInterfaceGeneric::getAllowCrossFlow() const { return well_ecl_.getAllowCrossFlow(); } -const Well& WellInterfaceGeneric::wellEcl() const +template +const Well& WellInterfaceGeneric::wellEcl() const { return well_ecl_; } -Well& WellInterfaceGeneric::wellEcl() +template +Well& WellInterfaceGeneric::wellEcl() { return well_ecl_; } -const PhaseUsage& WellInterfaceGeneric::phaseUsage() const +template +const PhaseUsage& WellInterfaceGeneric::phaseUsage() const { assert(phase_usage_ != nullptr); return *phase_usage_; } -double WellInterfaceGeneric::wsolvent() const +template +Scalar WellInterfaceGeneric::wsolvent() const { return wsolvent_; } -double WellInterfaceGeneric::rsRvInj() const +template +Scalar WellInterfaceGeneric::rsRvInj() const { return well_ecl_.getInjectionProperties().rsRvInj; } -void WellInterfaceGeneric::initInjMult(const std::vector& max_inj_mult) +template +void WellInterfaceGeneric:: +initInjMult(const std::vector& 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(max_inj_mult.size(), 1.); + this->inj_multiplier_ = std::vector(max_inj_mult.size(), 1.); } -void WellInterfaceGeneric::updateInjMult(std::vector& inj_multipliers, DeferredLogger& deferred_logger) const +template +void WellInterfaceGeneric:: +updateInjMult(std::vector& 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& inj_multipliers, D inj_multipliers = this->inj_multiplier_; } - - -double WellInterfaceGeneric::getInjMult(const int perf, - const double bhp, - const double perf_pres) const +template +Scalar WellInterfaceGeneric:: +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 +bool WellInterfaceGeneric:: +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& ws, - const double& simulationTime, - const bool& writeMessageToOPMLog, - WellTestState& wellTestState, - DeferredLogger& deferred_logger) const +template +void WellInterfaceGeneric:: +updateWellTestState(const SingleWellState& 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& ws // TODO: well can be shut/closed due to other reasons } -double WellInterfaceGeneric::getTHPConstraint(const SummaryState& summaryState) const +template +Scalar WellInterfaceGeneric:: +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 +bool WellInterfaceGeneric::underPredictionMode() const { return well_ecl_.predictionMode(); } -void WellInterfaceGeneric::initCompletions() +template +void WellInterfaceGeneric::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 +void WellInterfaceGeneric:: +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 +void WellInterfaceGeneric:: +setVFPProperties(const VFPProperties* vfp_properties_arg) { vfp_properties_ = vfp_properties_arg; } -void WellInterfaceGeneric::setPrevSurfaceRates(WellState& well_state, - const WellState& 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 +void WellInterfaceGeneric:: +setPrevSurfaceRates(WellState& well_state, + const WellState& 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 +void WellInterfaceGeneric:: +setGuideRate(const GuideRate* guide_rate_arg) { guide_rate_ = guide_rate_arg; } -void WellInterfaceGeneric::setWellEfficiencyFactor(const double efficiency_factor) +template +void WellInterfaceGeneric:: +setWellEfficiencyFactor(const Scalar efficiency_factor) { well_efficiency_factor_ = efficiency_factor; } -void WellInterfaceGeneric::setRepRadiusPerfLength() +template +void WellInterfaceGeneric::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 +void WellInterfaceGeneric:: +setWsolvent(const Scalar wsolvent) { wsolvent_ = wsolvent; } -void WellInterfaceGeneric::setDynamicThpLimit(const double thp_limit) +template +void WellInterfaceGeneric:: +setDynamicThpLimit(const Scalar thp_limit) { dynamic_thp_limit_ = thp_limit; } -std::optional WellInterfaceGeneric::getDynamicThpLimit() const +template +std::optional +WellInterfaceGeneric::getDynamicThpLimit() const { return dynamic_thp_limit_; } -void WellInterfaceGeneric::updatePerforatedCell(std::vector& is_cell_perforated) +template +void WellInterfaceGeneric:: +updatePerforatedCell(std::vector& is_cell_perforated) { - - for (int perf_idx = 0; perf_idx +bool WellInterfaceGeneric:: +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 +bool WellInterfaceGeneric::isOperableAndSolvable() const { return operability_status_.isOperableAndSolvable(); } -bool WellInterfaceGeneric::useVfpExplicit() const +template +bool WellInterfaceGeneric::useVfpExplicit() const { const auto& wvfpexp = well_ecl_.getWVFPEXP(); return (wvfpexp.explicit_lookup() || operability_status_.use_vfpexplicit); } -bool WellInterfaceGeneric::thpLimitViolatedButNotSwitched() const +template +bool WellInterfaceGeneric::thpLimitViolatedButNotSwitched() const { return operability_status_.thp_limit_violated_but_not_switched; } -double WellInterfaceGeneric::getALQ(const WellState& well_state) const +template +Scalar WellInterfaceGeneric:: +getALQ(const WellState& well_state) const { // no alq for injectors. if (isInjector()) @@ -514,8 +565,10 @@ double WellInterfaceGeneric::getALQ(const WellState& well_state) const return well_state.getALQ(name()); } -void WellInterfaceGeneric::reportWellSwitching(const SingleWellState &ws, - DeferredLogger& deferred_logger) const +template +void WellInterfaceGeneric:: +reportWellSwitching(const SingleWellState &ws, + DeferredLogger& deferred_logger) const { if (well_control_log_.empty()) return; @@ -534,7 +587,9 @@ void WellInterfaceGeneric::reportWellSwitching(const SingleWellState &ws } } -bool WellInterfaceGeneric::isPressureControlled(const WellState& well_state) const +template +bool WellInterfaceGeneric:: +isPressureControlled(const WellState& 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& well_st } } - - -bool WellInterfaceGeneric::wellUnderZeroRateTarget(const SummaryState& summary_state, - const WellState& well_state) const +template +bool WellInterfaceGeneric:: +wellUnderZeroRateTarget(const SummaryState& summary_state, + const WellState& 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& well_state) const +template +bool WellInterfaceGeneric:: +stopppedOrZeroRateTarget(const SummaryState& summary_state, + const WellState& well_state) const { return (this->wellIsStopped() || this->wellUnderZeroRateTarget(summary_state, well_state)); } -void WellInterfaceGeneric::resetWellOperability() +template +void WellInterfaceGeneric::resetWellOperability() { this->operability_status_.resetOperability(); } -double WellInterfaceGeneric::wmicrobes_() const +template +Scalar WellInterfaceGeneric::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 +Scalar WellInterfaceGeneric::wfoam_() const { auto injectorType = this->well_ecl_.injectorType(); @@ -603,7 +663,8 @@ double WellInterfaceGeneric::wfoam_() const } } -double WellInterfaceGeneric::wsalt_() const +template +Scalar WellInterfaceGeneric::wsalt_() const { auto injectorType = this->well_ecl_.injectorType(); @@ -616,13 +677,14 @@ double WellInterfaceGeneric::wsalt_() const } } -double WellInterfaceGeneric::woxygen_() const +template +Scalar WellInterfaceGeneric::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 +Scalar WellInterfaceGeneric::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 +Scalar WellInterfaceGeneric::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 +int WellInterfaceGeneric::polymerTable_() const { return this->well_ecl_.getPolymerProperties().m_skprpolytable; } -int WellInterfaceGeneric::polymerWaterTable_() const +template +int WellInterfaceGeneric::polymerWaterTable_() const { return this->well_ecl_.getPolymerProperties().m_skprwattable; } -int WellInterfaceGeneric::polymerInjTable_() const +template +int WellInterfaceGeneric::polymerInjTable_() const { return this->well_ecl_.getPolymerProperties().m_plymwinjtable; } -std::pair WellInterfaceGeneric:: -computeWellPotentials(std::vector& well_potentials, - const WellState& well_state) +template +std::pair WellInterfaceGeneric:: +computeWellPotentials(std::vector& well_potentials, + const WellState& well_state) { const int np = this->number_of_phases_; well_potentials.resize(np, 0.0); @@ -711,8 +779,8 @@ computeWellPotentials(std::vector& 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& well_potentials, return {compute_potential, bhp_controlled_well}; } -void WellInterfaceGeneric:: -checkNegativeWellPotentials(std::vector& well_potentials, +template +void WellInterfaceGeneric:: +checkNegativeWellPotentials(std::vector& 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& well_potentials, } } -void WellInterfaceGeneric:: +template +void WellInterfaceGeneric:: prepareForPotentialCalculations(const SummaryState& summary_state, - WellState& well_state, + WellState& well_state, Well::InjectionControls& inj_controls, Well::ProductionControls& prod_controls) const { @@ -793,4 +863,6 @@ prepareForPotentialCalculations(const SummaryState& summary_state, } } +template class WellInterfaceGeneric; + } // namespace Opm diff --git a/opm/simulators/wells/WellInterfaceGeneric.hpp b/opm/simulators/wells/WellInterfaceGeneric.hpp index fd0c10064..fa1d1eccc 100644 --- a/opm/simulators/wells/WellInterfaceGeneric.hpp +++ b/opm/simulators/wells/WellInterfaceGeneric.hpp @@ -47,6 +47,7 @@ template class SingleWellState; class Group; class Schedule; +template 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& rates) const; + void adaptRatesForVFP(std::vector& 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& well_state, - const WellState& prev_well_state) const; + void setPrevSurfaceRates(WellState& well_state, + const WellState& 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 getDynamicThpLimit() const; + void setWsolvent(const Scalar wsolvent); + void setDynamicThpLimit(const Scalar thp_limit); + std::optional getDynamicThpLimit() const; void updatePerforatedCell(std::vector& 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& perfDepth() const { return perf_depth_; } - const std::vector& perfDepth() const { - return perf_depth_; - } + std::vector& perfDepth() { return perf_depth_; } - std::vector& perfDepth() { - return perf_depth_; - } + const std::vector& wellIndex() const { return well_index_; } - const std::vector& wellIndex() const { - return well_index_; - } + const std::map>& getCompletions() const { return completions_; } - const std::map>& getCompletions() const { - return completions_; - } - - double getTHPConstraint(const SummaryState& summaryState) const; - double getALQ(const WellState& well_state) const; - double wsolvent() const; - double rsRvInj() const; + Scalar getTHPConstraint(const SummaryState& summaryState) const; + Scalar getALQ(const WellState& 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& max_inj_mult); + void initInjMult(const std::vector& 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& inj_multipliers, DeferredLogger& deferred_logger) const; + void updateInjMult(std::vector& 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& ws, DeferredLogger& deferred_logger) const; + void reportWellSwitching(const SingleWellState& 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& ws, + void updateWellTestState(const SingleWellState& ws, const double& simulationTime, const bool& writeMessageToOPMLog, WellTestState& wellTestState, DeferredLogger& deferred_logger) const; - bool isPressureControlled(const WellState& well_state) const; + bool isPressureControlled(const WellState& well_state) const; bool stopppedOrZeroRateTarget(const SummaryState& summary_state, - const WellState& well_state) const; + const WellState& 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& inj_fc_multiplier) + void updateFilterCakeMultipliers(const std::vector& 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& well_state) const; + const WellState& well_state) const; std::pair - computeWellPotentials(std::vector& well_potentials, - const WellState& well_state); + computeWellPotentials(std::vector& well_potentials, + const WellState& well_state); - void checkNegativeWellPotentials(std::vector& well_potentials, + void checkNegativeWellPotentials(std::vector& well_potentials, const bool checkOperability, DeferredLogger& deferred_logger); void prepareForPotentialCalculations(const SummaryState& summary_state, - WellState& well_state, + WellState& 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 ipr_a_; - mutable std::vector ipr_b_; + mutable std::vector ipr_a_; + mutable std::vector ipr_b_; // cell index for each well perforation std::vector well_cells_; // well index for each perforation - std::vector well_index_; + std::vector well_index_; // number of the perforations for this well int number_of_perforations_; // depth for each perforation - std::vector perf_depth_; + std::vector perf_depth_; // representative radius of the perforations, used in shear calculation - std::vector perf_rep_radius_; + std::vector perf_rep_radius_; // length of the perforations, use in shear calculation - std::vector perf_length_; + std::vector perf_length_; // well bore diameter - std::vector bore_diameters_; + std::vector bore_diameters_; /* * completions_ contains the mapping from completion id to connection indices @@ -364,7 +334,7 @@ protected: std::map> completions_; // reference depth for the BHP - double ref_depth_; + Scalar ref_depth_; // saturation table nubmer for each well perforation std::vector saturation_table_number_; @@ -373,25 +343,25 @@ protected: const PhaseUsage* phase_usage_; - double gravity_; - double wsolvent_; - std::optional dynamic_thp_limit_; + Scalar gravity_; + Scalar wsolvent_; + std::optional dynamic_thp_limit_; // recording the multiplier calculate from the keyword WINJMULT during the time step - mutable std::vector inj_multiplier_; + mutable std::vector 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 prev_inj_multiplier_; + std::vector prev_inj_multiplier_; // the multiplier due to injection filtration cake - std::vector inj_fc_multiplier_; + std::vector 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 well_control_log_; bool changed_to_open_this_step_ = true; }; diff --git a/opm/simulators/wells/WellInterfaceIndices.cpp b/opm/simulators/wells/WellInterfaceIndices.cpp index 2c34ddf5d..b0d672922 100644 --- a/opm/simulators/wells/WellInterfaceIndices.cpp +++ b/opm/simulators/wells/WellInterfaceIndices.cpp @@ -91,7 +91,7 @@ modelCompIdxToFlowCompIdx(const unsigned compIdx) const } template -double +typename WellInterfaceIndices::Scalar WellInterfaceIndices:: scalingFactor(const int phaseIdx) const { diff --git a/opm/simulators/wells/WellInterfaceIndices.hpp b/opm/simulators/wells/WellInterfaceIndices.hpp index 232c559f0..0671ba4b7 100644 --- a/opm/simulators/wells/WellInterfaceIndices.hpp +++ b/opm/simulators/wells/WellInterfaceIndices.hpp @@ -27,8 +27,7 @@ #include -namespace Opm -{ +namespace Opm { template class WellInterfaceIndices : public WellInterfaceFluidSystem @@ -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 Eval restrictEval(const EvalWell& in) const diff --git a/opm/simulators/wells/WellInterface_impl.hpp b/opm/simulators/wells/WellInterface_impl.hpp index ef8bdd627..b8011d0eb 100644 --- a/opm/simulators/wells/WellInterface_impl.hpp +++ b/opm/simulators/wells/WellInterface_impl.hpp @@ -89,10 +89,10 @@ namespace Opm void WellInterface:: init(const PhaseUsage* phase_usage_arg, - const std::vector& /* depth_arg */, - const double gravity_arg, + const std::vector& /* depth_arg */, + const Scalar gravity_arg, const int /* num_cells */, - const std::vector< Scalar >& B_avg, + const std::vector& B_avg, const bool changed_to_open_this_step) { this->phase_usage_ = phase_usage_arg; @@ -105,7 +105,7 @@ namespace Opm template - double + typename WellInterface::Scalar WellInterface:: wpolymer() const { @@ -121,7 +121,7 @@ namespace Opm template - double + typename WellInterface::Scalar WellInterface:: wfoam() const { @@ -135,7 +135,7 @@ namespace Opm template - double + typename WellInterface::Scalar WellInterface:: wsalt() const { @@ -147,7 +147,7 @@ namespace Opm } template - double + typename WellInterface::Scalar WellInterface:: wmicrobes() const { @@ -159,7 +159,7 @@ namespace Opm } template - double + typename WellInterface::Scalar WellInterface:: woxygen() const { @@ -177,7 +177,7 @@ namespace Opm // vtkblackoilmicpmodule.hh and ecloutputblackoilmodel.hh respectively). template - double + typename WellInterface::Scalar WellInterface:: wurea() const { @@ -273,7 +273,7 @@ namespace Opm const GroupState& 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 rates(this->num_components_); + std::vector 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(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(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 potentials; + std::vector 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(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(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 - std::optional + std::optional::Scalar> WellInterface:: 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:: solveWellWithBhp(const Simulator& simulator, const double dt, - const double bhp, + const Scalar bhp, WellState& well_state, DeferredLogger& deferred_logger) { @@ -750,9 +773,7 @@ namespace Opm const GroupState& 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 potentials; + std::vector 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 WellInterface::Scalar - WellInterface::volumetricSurfaceRateForConnection(int cellIdx, int phaseIdx) const { + WellInterface::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 bool WellInterface:: - gliftBeginTimeStepWellTestIterateWellEquations( - const Simulator& simulator, - const double dt, - WellState& well_state, - const GroupState& group_state, - DeferredLogger& deferred_logger) + gliftBeginTimeStepWellTestIterateWellEquations(const Simulator& simulator, + const double dt, + WellState& well_state, + const GroupState& 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 convert_coeff(this->number_of_phases_, 1.0); + std::vector 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; pcurrentStep()); - const double efficiencyFactor = well.getEfficiencyFactor(); - std::optional target = + const Scalar efficiencyFactor = well.getEfficiencyFactor(); + std::optional 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 fractions = initialWellRateFractions(simulator, well_state); + const std::vector fractions = initialWellRateFractions(simulator, well_state); double control_fraction = fractions[pu.phase_pos[Oil]]; if (control_fraction != 0.0) { for (int p = 0; p 0.0) { @@ -1213,11 +1236,11 @@ namespace Opm ws.surface_rates[p] *= controls.water_rate/current_rate; } } else { - const std::vector fractions = initialWellRateFractions(simulator, well_state); - double control_fraction = fractions[pu.phase_pos[Water]]; + const std::vector fractions = initialWellRateFractions(simulator, well_state); + const Scalar control_fraction = fractions[pu.phase_pos[Water]]; if (control_fraction != 0.0) { for (int p = 0; p 0.0) { @@ -1233,11 +1256,11 @@ namespace Opm ws.surface_rates[p] *= controls.gas_rate/current_rate; } } else { - const std::vector fractions = initialWellRateFractions(simulator, well_state); - double control_fraction = fractions[pu.phase_pos[Gas]]; + const std::vector fractions = initialWellRateFractions(simulator, well_state); + const Scalar control_fraction = fractions[pu.phase_pos[Gas]]; if (control_fraction != 0.0) { for (int p = 0; p 0.0) { @@ -1256,8 +1279,8 @@ namespace Opm ws.surface_rates[p] *= controls.liquid_rate/current_rate; } } else { - const std::vector fractions = initialWellRateFractions(simulator, well_state); - double control_fraction = fractions[pu.phase_pos[Water]] + fractions[pu.phase_pos[Oil]]; + const std::vector 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 convert_coeff(this->number_of_phases_, 1.0); + std::vector 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 fractions = initialWellRateFractions(simulator, well_state); + const std::vector fractions = initialWellRateFractions(simulator, well_state); for (int p = 0; p hrates(this->number_of_phases_,0.); + std::vector 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 hrates_resv(this->number_of_phases_,0.); + std::vector 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 fractions = initialWellRateFractions(simulator, well_state); + const std::vector fractions = initialWellRateFractions(simulator, well_state); for (int p = 0; padaptRatesForVFP(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 - std::vector + std::vector::Scalar> WellInterface:: initialWellRateFractions(const Simulator& simulator, const WellState& well_state) const { const int np = this->number_of_phases_; - std::vector scaling_factor(np); + std::vector 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; pnumber_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 well_q_s = computeCurrentWellRates(simulator, deferred_logger); + std::vector 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 - std::vector + std::vector::Scalar> WellInterface:: wellIndex(const int perf, const IntensiveQuantities& intQuants, - const double trans_mult, - const SingleWellState& ws) const + const Scalar trans_mult, + const SingleWellState& 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:: updateConnectionDFactor(const Simulator& simulator, - SingleWellState& ws) const + SingleWellState& ws) const { if (! this->well_ecl_.getWDFAC().useDFactor()) { return; @@ -1600,11 +1621,11 @@ namespace Opm } template - double + typename WellInterface::Scalar WellInterface:: computeConnectionDFactor(const int perf, const IntensiveQuantities& intQuants, - const SingleWellState& ws) const + const SingleWellState& 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:: updateConnectionTransmissibilityFactor(const Simulator& simulator, - SingleWellState& ws) const + SingleWellState& 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 rates(this->number_of_phases_, 0.0); + std::vector 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:: computeConnLevelProdInd(const FluidState& fs, - const std::function& connPICalc, + const std::function& connPICalc, const std::vector& 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:: computeConnLevelInjInd(const FluidState& fs, const Phase preferred_phase, - const std::function& connIICalc, + const std::function& connIICalc, const std::vector& mobility, - double* connII, + Scalar* connII, DeferredLogger& deferred_logger) const { // Assumes single phase injection diff --git a/opm/simulators/wells/WellTest.hpp b/opm/simulators/wells/WellTest.hpp index c16ab3b06..d255734e3 100644 --- a/opm/simulators/wells/WellTest.hpp +++ b/opm/simulators/wells/WellTest.hpp @@ -34,14 +34,14 @@ class DeferredLogger; struct PhaseUsage; template class SingleWellState; class WellEconProductionLimits; -class WellInterfaceGeneric; +template 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& well) : well_(well) {} void updateWellTestStateEconomic(const SingleWellState& ws, const double simulation_time, @@ -95,7 +95,7 @@ private: DeferredLogger& deferred_logger) const; - const WellInterfaceGeneric& well_; //!< Reference to well interface + const WellInterfaceGeneric& well_; //!< Reference to well interface }; }