WellInterfaceGeneric.hpp

Go to the documentation of this file.

/*
  Copyright 2017 SINTEF Digital, Mathematics and Cybernetics.
  Copyright 2017 Statoil ASA.
  Copyright 2017 IRIS
  Copyright 2019 Norce
 
  This file is part of the Open Porous Media project (OPM).
 
  OPM is free software: you can redistribute it and/or modify
  it under the terms of the GNU General Public License as published by
  the Free Software Foundation, either version 3 of the License, or
  (at your option) any later version.
 
  OPM is distributed in the hope that it will be useful,
  but WITHOUT ANY WARRANTY; without even the implied warranty of
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  GNU General Public License for more details.
 
  You should have received a copy of the GNU General Public License
  along with OPM.  If not, see <http://www.gnu.org/licenses/>.
*/
 
 
#ifndef OPM_WELLINTERFACE_GENERIC_HEADER_INCLUDED
#define OPM_WELLINTERFACE_GENERIC_HEADER_INCLUDED
 
#include <dune/common/timer.hh>
 
#include <opm/input/eclipse/Schedule/Well/Well.hpp>
#include <opm/simulators/flow/BlackoilModelParameters.hpp>
#include <opm/simulators/wells/RuntimePerforation.hpp>
#include <opm/simulators/wells/WellIndexFracture.hpp>
 
#include <ctime>
#include <map>
#include <memory>
#include <optional>
#include <string>
#include <vector>
 
namespace Opm
{
 
class DeferredLogger;
class GuideRate;
class UnitSystem;
template<class Scalar> class ParallelWellInfo;
template<class Scalar> struct PerforationData;
class SummaryState;
template<typename Scalar, typename IndexTraits> class VFPProperties;
class WellTestState;
template<typename Scalar, typename IndexTraits> class WellState;
template<typename Scalar, typename IndexTraits> class SingleWellState;
class Group;
class Schedule;
template<typename IndexTraits> class PhaseUsageInfo;
 
template<typename Scalar, typename IndexTraits>
class WellInterfaceGeneric {
public:
    using ModelParameters = BlackoilModelParameters<Scalar>;
 
    using WellStateType = WellState<Scalar, IndexTraits>;
 
    WellInterfaceGeneric(const Well& well,
                         const ParallelWellInfo<Scalar>& parallel_well_info,
                         const int time_step,
                         const ModelParameters& param,
                         const int pvtRegionIdx,
                         const int num_conservation_quantities,
                         const int num_phases,
                         const int index_of_well,
                         const PhaseUsageInfo<IndexTraits>& phase_usage,
                         const std::vector<PerforationData<Scalar>>& perf_data);
 
    virtual ~WellInterfaceGeneric() = default;
 
    const std::vector<PerforationData<Scalar>>& perforationData() const;
 
    const std::string& name() const;
 
    bool isInjector() const;
 
    bool isProducer() const;
 
    const std::vector<int>& cells() const { return well_cells_; }
 
    int indexOfWell() const;
 
    void adaptRatesForVFP(std::vector<Scalar>& rates) const;
 
    const Well& wellEcl() const;
    Well& wellEcl();
    const PhaseUsageInfo<IndexTraits>& phaseUsage() const;
 
    bool underPredictionMode() const;
 
    // whether the well is operable
    bool isOperableAndSolvable() const;
    bool useVfpExplicit () const;
    bool thpLimitViolatedButNotSwitched() const;
 
    void initCompletions();
    void closeCompletions(const WellTestState& wellTestState);
 
    void setVFPProperties(const VFPProperties<Scalar, IndexTraits>* vfp_properties_arg);
    void setPrevSurfaceRates(WellStateType& well_state,
                             const WellStateType& prev_well_state) const;
    void setGuideRate(const GuideRate* guide_rate_arg);
    void setWellEfficiencyFactor(const Scalar efficiency_factor);
    void setRepRadiusPerfLength();
    void setWsolvent(const Scalar wsolvent);
    void setDynamicThpLimit(const Scalar thp_limit);
    std::optional<Scalar> getDynamicThpLimit() const;
    void setDynamicThpLimit(const std::optional<Scalar> thp_limit);
    void updatePerforatedCell(std::vector<bool>& is_cell_perforated);
 
    bool wellHasTHPConstraints(const SummaryState& summaryState) const;
 
    void stopWell() { this->wellStatus_ = Well::Status::STOP; }
    void openWell() { this->wellStatus_ = Well::Status::OPEN; }
    Well::Status wellStatus() { return this->wellStatus_;}
 
    bool wellIsStopped() const { return this->wellStatus_ == Well::Status::STOP; }
 
    int currentStep() const { return this->current_step_; }
 
    int pvtRegionIdx() const { return pvtRegionIdx_; }
 
    const GuideRate* guideRate() const { return guide_rate_; }
 
    int numConservationQuantities() const { return num_conservation_quantities_; }
 
    int numPhases() const { return number_of_phases_; }
 
    int numLocalPerfs() const { return number_of_local_perforations_; }
 
    Scalar refDepth() const { return ref_depth_; }
 
    Scalar gravity() const { return gravity_; }
 
    const VFPProperties<Scalar, IndexTraits>* vfpProperties() const { return vfp_properties_; }
 
    const ParallelWellInfo<Scalar>& parallelWellInfo() const { return parallel_well_info_; }
 
    const std::vector<Scalar>& perfDepth() const { return perf_depth_; }
 
    std::vector<Scalar>& perfDepth() { return perf_depth_; }
 
    const std::vector<Scalar>& wellIndex() const { return well_index_; }
 
    std::vector<Scalar>& wellIndex() { return well_index_; }
 
    const std::map<int,std::vector<int>>& getCompletions() const { return completions_; }
 
    Scalar getTHPConstraint(const SummaryState& summaryState) const;
    Scalar getALQ(const WellStateType& 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<Scalar>& max_inj_mult);
 
    // update the InjMult information at the end of the time step, so it can be used for later.
    void updateInjMult(std::vector<Scalar>& inj_multipliers,
                       DeferredLogger& deferred_logger) const;
 
    // Note:: for multisegment wells, bhp is actually segment pressure in practice based on observation
    // it might change in the future
    Scalar getInjMult(const int local_perf_index, const Scalar bhp, const Scalar perf_pres, DeferredLogger& dlogger) 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<Scalar, IndexTraits>& ws,
                             DeferredLogger& deferred_logger) const;
 
    bool changedToOpenThisStep() const { return this->changed_to_open_this_step_; }
 
    void updateWellTestState(const SingleWellState<Scalar, IndexTraits>& ws,
                             const double& simulationTime,
                             const bool& writeMessageToOPMLog,
                             const bool during_well_test,
                             const bool zero_group_target,
                             WellTestState& wellTestState,
                             const UnitSystem& unit_system,
                             const std::time_t start_time,
                             DeferredLogger& deferred_logger,
                             std::string* closure_reason = nullptr) const;
 
    bool isPressureControlled(const WellStateType& well_state) const;
 
    static void applyWeldrawRateLimit(const Well& well_ecl,
                                      const std::optional<Scalar>& weldraw_max_rate,
                                      Well::ProductionControls& controls);
 
    Well::ProductionControls
    productionControlsWithWeldraw(const SummaryState& summary_state,
                                  const SingleWellState<Scalar, IndexTraits>& ws) const;
 
    Scalar wellEfficiencyFactor() const { return well_efficiency_factor_; }
 
    void updateFilterCakeMultipliers(const std::vector<Scalar>& inj_fc_multiplier)
    {
        inj_fc_multiplier_ = inj_fc_multiplier;
    }
 
    void resetWellOperability();
 
    virtual std::vector<Scalar> getPrimaryVars() const
    {
        return {};
    }
 
    virtual int setPrimaryVars(typename std::vector<Scalar>::const_iterator)
    {
        return 0;
    }
 
    virtual Scalar connectionDensity(const int globalConnIdx,
                                     const int openConnIdx) const = 0;
 
    virtual Scalar totalReservoirVoidageRate(const std::vector<Scalar>& surface_rates) const = 0;
 
    void addPerforations(const std::vector<RuntimePerforation>& perfs);
 
    void addFracturePerforations(const std::vector<RuntimePerforation>& perfs);
 
    const std::vector<WellIndexFracture>& wellIndexFracture() const
    { return well_index_fracture_; }
 
    struct WellSolveStats
    {
        double solve_time{0.0};            
        double potential_solve_time{0.0};  
        double assemble_time{0.0};         
        double linear_solve_time{0.0};     
        unsigned int iterations{0};            
        unsigned int potential_iterations{0};  
    };
 
    const WellSolveStats& solveStats() const { return *solve_stats_; }
    void resetSolveStats() const { *solve_stats_ = WellSolveStats{}; }
 
protected:
    class WellSolveScope
    {
    public:
        WellSolveScope(const WellInterfaceGeneric& well, const int& iterations)
            : well_(well), iterations_(iterations)
        {
            well_.in_well_solve_ = true;
            timer_.start();
        }
        ~WellSolveScope()
        {
            well_.in_well_solve_ = false;
            auto& stats = *well_.solve_stats_;
            if (well_.computing_well_potentials_) {
                stats.potential_solve_time += timer_.stop();
                stats.potential_iterations += iterations_;
            } else {
                stats.solve_time += timer_.stop();
                stats.iterations += iterations_;
            }
        }
    private:
        const WellInterfaceGeneric& well_;
        const int& iterations_;
        Dune::Timer timer_{};
    };
 
    class SolvePartTimer
    {
    public:
        SolvePartTimer(const WellInterfaceGeneric& well, double& target)
            : well_(well), target_(target)
        {
            timer_.start();
        }
        ~SolvePartTimer()
        {
            if (well_.in_well_solve_) {
                target_ += timer_.stop();
            }
        }
    private:
        const WellInterfaceGeneric& well_;
        double& target_;
        Dune::Timer timer_{};
    };
 
    class PotentialCalculationScope
    {
    public:
        explicit PotentialCalculationScope(const WellInterfaceGeneric& well)
            : well_(well), previous_(well.computing_well_potentials_)
        {
            well_.computing_well_potentials_ = true;
        }
        ~PotentialCalculationScope()
        {
            well_.computing_well_potentials_ = previous_;
        }
    private:
        const WellInterfaceGeneric& well_;
        bool previous_;
    };
 
    WellSolveScope solveScope(const int& iterations) const
    { return WellSolveScope(*this, iterations); }
 
    SolvePartTimer solveAssembleTimer() const
    { return SolvePartTimer(*this, solve_stats_->assemble_time); }
 
    SolvePartTimer solveLinearSolveTimer() const
    { return SolvePartTimer(*this, solve_stats_->linear_solve_time); }
 
    PotentialCalculationScope potentialCalculationScope() const
    { return PotentialCalculationScope(*this); }
 
    bool getAllowCrossFlow() 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 wellUnderZeroRateTargetIndividual(const SummaryState& summary_state,
                                           const WellState<Scalar, IndexTraits>& well_state) const;
 
    bool wellUnderGroupControl(const SingleWellState<Scalar, IndexTraits>& ws) const;
 
    std::pair<bool,bool>
    computeWellPotentials(std::vector<Scalar>& well_potentials,
                          const WellStateType& well_state);
 
    void checkNegativeWellPotentials(std::vector<Scalar>& well_potentials,
                                     const bool checkOperability,
                                     DeferredLogger& deferred_logger);
 
    // Remove all other controls than THP and BHP
    void onlyKeepBHPandTHPcontrols(const SummaryState& summary_state,
                                   WellStateType& well_state,
                                   Well::InjectionControls& inj_controls,
                                   Well::ProductionControls& prod_controls) const;
 
    void resetDampening()
    {
        std::ranges::fill(this->inj_multiplier_damp_factor_, 1.0);
    }
 
    void updateGroupTargetFallbackFlag(WellState<Scalar, IndexTraits>& well_state,
                                       const std::vector<Scalar>& scaled_well_fractions,
                                       DeferredLogger& deferred_logger) const;
 
    // definition of the struct OperabilityStatus
    struct OperabilityStatus
    {
        bool isOperableAndSolvable() const
        {
            if (!operable_under_only_bhp_limit || !solvable || has_negative_potentials) {
                return false;
            }
 
            return isOperableUnderBHPLimit() || isOperableUnderTHPLimit();
        }
 
        bool isOperableUnderBHPLimit() const
        {
            return operable_under_only_bhp_limit && obey_thp_limit_under_bhp_limit;
        }
 
        bool isOperableUnderTHPLimit() const
        {
            return can_obtain_bhp_with_thp_limit && obey_bhp_limit_with_thp_limit;
        }
 
        void resetOperability()
        {
            operable_under_only_bhp_limit = true;
            obey_thp_limit_under_bhp_limit = true;
            can_obtain_bhp_with_thp_limit = true;
            obey_bhp_limit_with_thp_limit = true;
        }
 
        // whether the well can be operated under bhp limit
        // without considering other limits.
        // if it is false, then the well is not operable for sure.
        bool operable_under_only_bhp_limit = true;
        // if the well can be operated under bhp limit, will it obey(not violate)
        // the thp limit when operated under bhp limit
        bool obey_thp_limit_under_bhp_limit = true;
        // whether the well operate under the thp limit only
        bool can_obtain_bhp_with_thp_limit = true;
        // whether the well obey bhp limit when operated under thp limit
        bool obey_bhp_limit_with_thp_limit = true;
        // the well is solveable
        bool solvable = true;
        // the well have non positive potentials
        bool has_negative_potentials = false;
        //thp limit violated but not switched
        mutable bool thp_limit_violated_but_not_switched = false;
 
        bool use_vfpexplicit = false;
    };
 
    OperabilityStatus operability_status_;
 
    Well well_ecl_;
 
    const ParallelWellInfo<Scalar>& parallel_well_info_;
    const int current_step_;
    const ModelParameters& param_;
 
    // The pvt region of the well. We assume
    // We assume a well to not penetrate more than one pvt region.
    const int pvtRegionIdx_;
 
    // TODO: with energy conservation joins, we probably should change this to be num_mass_conservation_quantities_
    // or something like that. Energy term calculate somehow differently from the mass term
    const int num_conservation_quantities_;
 
    // number of phases
    int number_of_phases_;
 
    // the index of well in Wells struct
    int index_of_well_;
 
    const PhaseUsageInfo<IndexTraits>& phase_usage_;
 
    const std::vector<PerforationData<Scalar>>* perf_data_;
 
    // the vectors used to describe the inflow performance relationship (IPR)
    // 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<Scalar> ipr_a_;
    mutable std::vector<Scalar> ipr_b_;
 
    // cell index for each well perforation
    std::vector<int> well_cells_;
 
    // well index for each perforation
    std::vector<Scalar> well_index_;
 
    // fracture contribution to the well index per perforation; empty when no
    // fracture model has registered contributions (see addFracturePerforations)
    std::vector<WellIndexFracture> well_index_fracture_;
 
    // number of the perforations for this well on this process
    int number_of_local_perforations_;
 
    // depth for each perforation
    std::vector<Scalar> perf_depth_;
 
    // representative radius of the perforations, used in shear calculation
    std::vector<Scalar> perf_rep_radius_;
 
    // length of the perforations, use in shear calculation
    std::vector<Scalar> perf_length_;
 
    // well bore diameter
    std::vector<Scalar> bore_diameters_;
 
    /*
     *  completions_ contains the mapping from completion id to connection indices
     *  {
     *      2 : [ConnectionIndex, ConnectionIndex],
     *      1 : [ConnectionIndex, ConnectionIndex, ConnectionIndex],
     *      5 : [ConnectionIndex],
     *      7 : [ConnectionIndex]
     *      ...
     *   }
     *   The integer IDs correspond to the COMPLETION id given by the COMPLUMP keyword.
     *   When there is no COMPLUMP keyword used, a default completion number will be assigned
     *   based on the order of the declaration of the connections.
     *   Since the connections not OPEN is not included in the Wells, so they will not be considered
     *   in this mapping relation.
     */
    std::map<int, std::vector<int>> completions_;
 
    // reference depth for the BHP
    Scalar ref_depth_;
 
    // saturation table nubmer for each well perforation
    std::vector<int> saturation_table_number_;
 
    Well::Status wellStatus_;
 
    Scalar gravity_;
    Scalar wsolvent_;
    std::optional<Scalar> dynamic_thp_limit_;
 
    // recording the multiplier calculate from the keyword WINJMULT during the time step
    mutable std::vector<Scalar> inj_multiplier_;
 
    // the injection multiplier from the previous running, it is mostly used for CIRR mode
    // which intends to keep the fracturing open
    std::vector<Scalar> prev_inj_multiplier_;
 
    // WINJMULT multipliers for previous iteration (used for oscillation detection)
    mutable std::vector<Scalar> inj_multiplier_previter_;
    // WINJMULT dampening factors (used in case of oscillations)
    mutable std::vector<Scalar> inj_multiplier_damp_factor_;
 
    // the multiplier due to injection filtration cake
    std::vector<Scalar> inj_fc_multiplier_;
 
    Scalar well_efficiency_factor_;
    const VFPProperties<Scalar, IndexTraits>* vfp_properties_;
    const GuideRate* guide_rate_;
 
    std::vector<std::string> well_control_log_;
 
    bool changed_to_open_this_step_ = true;
 
    // statistics for standalone well solves; shared so that well copies used
    // during e.g. potential calculations contribute to the same statistics
    std::shared_ptr<WellSolveStats> solve_stats_ = std::make_shared<WellSolveStats>();
    // true while a well potential calculation is in progress
    mutable bool computing_well_potentials_ = false;
    // true while a standalone well solve is in progress
    mutable bool in_well_solve_ = false;
};
 
}
 
#endif // OPM_WELLINTERFACE_GENERIC_HEADER_INCLUDED