WellInterface.hpp

/*
  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_HEADER_INCLUDED
#define OPM_WELLINTERFACE_HEADER_INCLUDED

// NOTE: GasLiftSingleWell.hpp includes StandardWell.hpp which includes ourself
//   (WellInterface.hpp), so we need to forward declare GasLiftSingleWell
//   for it to be defined in this file. Similar for BlackoilWellModel
namespace Opm {
    template<typename TypeTag> class GasLiftSingleWell;
    template<typename TypeTag> class BlackoilWellModel;
}

#include <opm/common/OpmLog/OpmLog.hpp>
#include <opm/common/ErrorMacros.hpp>
#include <opm/common/Exceptions.hpp>

#include <opm/input/eclipse/Schedule/Well/WellTestState.hpp>

#include <opm/material/fluidstates/BlackOilFluidState.hpp>

#include <opm/models/blackoil/blackoilproperties.hh>

#include <opm/simulators/linalg/linalgproperties.hh>

#include <opm/simulators/wells/BlackoilWellModel.hpp>
#include <opm/simulators/wells/GasLiftGroupInfo.hpp>
#include <opm/simulators/wells/GasLiftSingleWell.hpp>
#include <opm/simulators/wells/GasLiftSingleWellGeneric.hpp>
#include <opm/simulators/wells/PerforationData.hpp>
#include <opm/simulators/wells/WellInterfaceIndices.hpp>
#include <opm/simulators/wells/WellProdIndexCalculator.hpp>
#include <opm/simulators/wells/WellState.hpp>

#include <opm/simulators/timestepping/ConvergenceReport.hpp>

#include <opm/simulators/utils/DeferredLogger.hpp>

#include <dune/common/fmatrix.hh>
#include <dune/istl/bcrsmatrix.hh>
#include <dune/istl/matrixmatrix.hh>

#include <opm/material/densead/Evaluation.hpp>

#include <limits>
#include <vector>

namespace Opm
{

class WellInjectionProperties;
class WellProductionProperties;
template<typename Scalar, typename IndexTraits> class GroupStateHelper;

template<typename TypeTag>
class WellInterface : public WellInterfaceIndices<GetPropType<TypeTag, Properties::FluidSystem>,
                                                  GetPropType<TypeTag, Properties::Indices>>
{
    using Base = WellInterfaceIndices<GetPropType<TypeTag, Properties::FluidSystem>,
                                      GetPropType<TypeTag, Properties::Indices>>;
public:
    using Grid = GetPropType<TypeTag, Properties::Grid>;
    using Simulator = GetPropType<TypeTag, Properties::Simulator>;
    using Scalar = GetPropType<TypeTag, Properties::Scalar>;
    using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
    using IndexTraits = typename FluidSystem::IndexTraitsType;
    using Indices = GetPropType<TypeTag, Properties::Indices>;
    using IntensiveQuantities = GetPropType<TypeTag, Properties::IntensiveQuantities>;
    using MaterialLaw = GetPropType<TypeTag, Properties::MaterialLaw>;
    using SparseMatrixAdapter = GetPropType<TypeTag, Properties::SparseMatrixAdapter>;
    using RateVector = GetPropType<TypeTag, Properties::RateVector>;
    using GasLiftSingleWell = ::Opm::GasLiftSingleWell<TypeTag>;
    using GLiftEclWells = typename GasLiftGroupInfo<Scalar, IndexTraits>::GLiftEclWells;

    using VectorBlockType = Dune::FieldVector<Scalar, Indices::numEq>;
    using MatrixBlockType = Dune::FieldMatrix<Scalar, Indices::numEq, Indices::numEq>;
    using Eval = typename Base::Eval;
    using BVector = Dune::BlockVector<VectorBlockType>;
    using PressureMatrix = Dune::BCRSMatrix<Opm::MatrixBlock<Scalar, 1, 1>>;

    using WellStateType = WellState<Scalar, IndexTraits>;
    using SingleWellStateType = SingleWellState<Scalar, IndexTraits>;
    using GroupStateHelperType = GroupStateHelper<Scalar, IndexTraits>;

    using RateConverterType =
    typename WellInterfaceFluidSystem<FluidSystem>::RateConverterType;

    using WellInterfaceFluidSystem<FluidSystem>::Gas;
    using WellInterfaceFluidSystem<FluidSystem>::Oil;
    using WellInterfaceFluidSystem<FluidSystem>::Water;

    using ModelParameters = typename Base::ModelParameters;

    static constexpr bool has_solvent = getPropValue<TypeTag, Properties::EnableSolvent>();
    static constexpr bool has_zFraction = getPropValue<TypeTag, Properties::EnableExtbo>();
    static constexpr bool has_polymer = getPropValue<TypeTag, Properties::EnablePolymer>();
    static constexpr EnergyModules energyModuleType = getPropValue<TypeTag, Properties::EnergyModuleType>();
    static constexpr bool has_energy = energyModuleType == EnergyModules::FullyImplicitThermal;

    // flag for polymer molecular weight related
    static constexpr bool has_polymermw = getPropValue<TypeTag, Properties::EnablePolymerMW>();
    static constexpr bool has_foam = getPropValue<TypeTag, Properties::EnableFoam>();
    static constexpr bool has_brine = getPropValue<TypeTag, Properties::EnableBrine>();
    static constexpr bool has_watVapor = getPropValue<TypeTag, Properties::EnableVapwat>();
    static constexpr bool has_disgas_in_water = getPropValue<TypeTag, Properties::EnableDisgasInWater>();
    static constexpr bool has_saltPrecip = getPropValue<TypeTag, Properties::EnableSaltPrecipitation>();
    static constexpr bool has_bioeffects = getPropValue<TypeTag, Properties::EnableBioeffects>();
    static constexpr bool has_micp = Indices::enableMICP;

    template<class ValueType>
    using BlackOilFluidStateType = BlackOilFluidState<ValueType,
                                                      FluidSystem,
                                                      energyModuleType != EnergyModules::NoTemperature,
                                                      energyModuleType == EnergyModules::FullyImplicitThermal,
                                                      Indices::compositionSwitchIdx != std::numeric_limits<unsigned>::max(),
                                                      has_watVapor,
                                                      has_brine,
                                                      has_saltPrecip,
                                                      has_disgas_in_water,
                                                      has_solvent,
                                                      Indices::numPhases>;

    // fluid state for the reservoir fluid
    using FluidState = BlackOilFluidStateType<Eval>;

    WellInterface(const Well& well,
                  const ParallelWellInfo<Scalar>& pw_info,
                  const int time_step,
                  const ModelParameters& param,
                  const RateConverterType& rate_converter,
                  const int pvtRegionIdx,
                  const int num_conservation_quantities,
                  const int num_phases,
                  const int index_of_well,
                  const std::vector<PerforationData<Scalar>>& perf_data);

    virtual ~WellInterface() = default;

    virtual void init(const std::vector<Scalar>& depth_arg,
                      const Scalar gravity_arg,
                      const std::vector<Scalar>& B_avg,
                      const bool changed_to_open_this_step);

    virtual ConvergenceReport getWellConvergence(const GroupStateHelperType& groupStateHelper,
                                                 const std::vector<Scalar>& B_avg,
                                                 const bool relax_tolerance) const = 0;

    virtual void solveEqAndUpdateWellState(const Simulator& simulator,
                                           const GroupStateHelperType& groupStateHelper,
                                           WellStateType& well_state) = 0;

    void assembleWellEq(const Simulator& simulator,
                        const double dt,
                        const GroupStateHelperType& groupStateHelper,
                        WellStateType& well_state);

    void assembleWellEqWithoutIteration(const Simulator& simulator,
                                        const GroupStateHelperType& groupStateHelper,
                                        const double dt,
                                        WellStateType& well_state,
                                        const bool solving_with_zero_rate);

    // Convenience overload that gets scaled fractions internally
    void updateGroupTargetFallbackFlag(WellStateType& well_state,
                                       DeferredLogger& deferred_logger) const;

    // TODO: better name or further refactoring the function to make it more clear
    void prepareWellBeforeAssembling(const Simulator& simulator,
                                     const double dt,
                                     const GroupStateHelperType& groupStateHelper,
                                     WellStateType& well_state);

    virtual void computeWellRatesWithBhp(const Simulator& ebosSimulator,
                                         const Scalar& bhp,
                                         std::vector<Scalar>& well_flux,
                                         DeferredLogger& deferred_logger) const = 0;

    virtual std::optional<Scalar>
    computeBhpAtThpLimitProdWithAlq(const Simulator& ebos_simulator,
                                    const GroupStateHelperType& groupStateHelper,
                                    const SummaryState& summary_state,
                                    const Scalar alq_value,
                                    bool iterate_if_no_solution) const = 0;

    std::optional<Scalar>
    computeBhpAtThpLimitProdWithAlqUsingIPR(const Simulator& simulator,
                                            const WellStateType& well_state,
                                            Scalar bhp,
                                            const SummaryState& summary_state,
                                            const Scalar alq_value);
    virtual void recoverWellSolutionAndUpdateWellState(const Simulator& simulator,
                                                       const BVector& x,
                                                       const GroupStateHelperType& groupStateHelper,
                                                       WellStateType& well_state) = 0;

    virtual void apply(const BVector& x, BVector& Ax) const = 0;

    virtual void apply(BVector& r) const = 0;

    // TODO: before we decide to put more information under mutable, this function is not const
    virtual void computeWellPotentials(const Simulator& simulator,
                                       const WellStateType& well_state,
                                       const GroupStateHelperType& groupStateHelper,
                                       std::vector<Scalar>& well_potentials) = 0;

    virtual void updateWellStateWithTarget(const Simulator& simulator,
                                           const GroupStateHelperType& groupStateHelper,
                                           WellStateType& well_state) const;

    virtual void scaleSegmentRatesAndPressure(WellStateType& well_state) const;

    virtual void computeWellRatesWithBhpIterations(const Simulator& simulator,
                                                   const Scalar& bhp,
                                                   const GroupStateHelperType& groupStateHelper,
                                                   std::vector<Scalar>& well_flux) const = 0;

    bool wellUnderZeroRateTarget(const GroupStateHelperType& groupStateHelper) const;

    bool stoppedOrZeroRateTarget(const GroupStateHelperType& groupStateHelper) const;

    bool updateWellStateWithTHPTargetProd(const Simulator& simulator,
                                          WellStateType& well_state,
                                          const GroupStateHelperType& groupStateHelper) const;

    bool wellUnderZeroGroupRateTarget(const GroupStateHelperType& groupStateHelper,
                                      const std::optional<bool> group_control = std::nullopt) const;

    enum class IndividualOrGroup { Individual, Group, Both };
    bool updateWellControl(const Simulator& simulator,
                           const IndividualOrGroup iog,
                           const GroupStateHelperType& groupStateHelper,
                           WellStateType& well_state) /* const */;

    bool updateWellControlAndStatusLocalIteration(const Simulator& simulator,
                                                  const GroupStateHelperType& groupStateHelper,
                                                  const Well::InjectionControls& inj_controls,
                                                  const Well::ProductionControls& prod_controls,
                                                  const Scalar WQTotal,
                                                  WellStateType& well_state,
                                                  const bool fixed_control,
                                                  const bool fixed_status,
                                                  const bool solving_with_zero_rate);

    virtual void updatePrimaryVariables(const GroupStateHelperType& groupStateHelper) = 0;

    virtual void calculateExplicitQuantities(const Simulator& simulator,
                                             const GroupStateHelperType& groupStateHelper) = 0; // should be const?

    virtual void updateProductivityIndex(const Simulator& simulator,
                                         const WellProdIndexCalculator<Scalar>& wellPICalc,
                                         WellStateType& well_state,
                                         DeferredLogger& deferred_logger) const = 0;

    // Add well contributions to matrix
    virtual void addWellContributions(SparseMatrixAdapter&) const = 0;

    virtual void addWellPressureEquations(PressureMatrix& mat,
                                          const BVector& x,
                                          const int pressureVarIndex,
                                          const bool use_well_weights,
                                          const WellStateType& well_state) const = 0;

    void addCellRates(std::map<int, RateVector>& cellRates_) const;

    Scalar volumetricSurfaceRateForConnection(int cellIdx, int phaseIdx) const;

    // TODO: theoretically, it should be a const function
    // Simulator is not const is because that assembleWellEq is non-const Simulator
    void wellTesting(const Simulator& simulator,
                     const double simulation_time,
                     const GroupStateHelperType& groupStateHelper,
                     WellStateType& well_state,
                     WellTestState& welltest_state,
                     GLiftEclWells& ecl_well_map,
                     std::map<std::string, double>& open_times);

    void checkWellOperability(const Simulator& simulator,
                              const WellStateType& well_state,
                              const GroupStateHelperType& groupStateHelper);

    void gliftBeginTimeStepWellTestUpdateALQ(const Simulator& simulator,
                                             WellStateType& well_state,
                                             const GroupState<Scalar>& group_state,
                                             GLiftEclWells& ecl_well_map,
                                             DeferredLogger& deferred_logger);

    // check whether the well is operable under the current reservoir condition
    // mostly related to BHP limit and THP limit
    void updateWellOperability(const Simulator& simulator,
                               const WellStateType& well_state,
                               const GroupStateHelperType& groupStateHelper);

    bool updateWellOperabilityFromWellEq(const Simulator& simulator,
                                         const GroupStateHelperType& groupStateHelper);

    // update perforation water throughput based on solved water rate
    virtual void updateWaterThroughput(const double dt,
                                       WellStateType& well_state) const = 0;

    virtual std::vector<Scalar>
    computeCurrentWellRates(const Simulator& simulator,
                            DeferredLogger& deferred_logger) const = 0;

    virtual void getScaledWellFractions(std::vector<Scalar>& scaled_fractions,
                                        DeferredLogger& deferred_logger) const = 0;

    void initializeProducerWellState(const Simulator& simulator,
                                     WellStateType& well_state,
                                     DeferredLogger& deferred_logger) const;

    void solveWellEquation(const Simulator& simulator,
                           const GroupStateHelperType& groupStateHelper,
                           WellStateType& well_state);

    const std::vector<RateVector>& connectionRates() const
    {
        return connectionRates_;
    }

    void updateConnectionDFactor(const Simulator& simulator,
                                 SingleWellStateType& ws) const;

    void updateConnectionTransmissibilityFactor(const Simulator& simulator,
                                                SingleWellStateType& ws) const;

    virtual bool iterateWellEqWithSwitching(const Simulator& simulator,
                                            const double dt,
                                            const WellInjectionControls& inj_controls,
                                            const WellProductionControls& prod_controls,
                                            const GroupStateHelperType& groupStateHelper,
                                            WellStateType& well_state,
                                            const bool fixed_control,
                                            const bool fixed_status,
                                            const bool solving_with_zero_rate) = 0;

    virtual Scalar maxPerfPress(const Simulator& simulator) const = 0;

    virtual void updateIPRImplicit(const Simulator& simulator,
                                   const GroupStateHelperType& groupStateHelper,
                                   WellStateType& well_state) = 0;

    static constexpr int numResDofs = Indices::numEq;
    static constexpr int numWellDofs = numResDofs + 1;  // NB will fail for for thermal for now
    using BMatrix = Dune::BCRSMatrix<Dune::FieldMatrix<Scalar, numWellDofs, numResDofs>>;
    using CMatrix = Dune::BCRSMatrix<Dune::FieldMatrix<Scalar, numResDofs, numWellDofs>>;
    using DMatrix = Dune::BCRSMatrix<Dune::FieldMatrix<Scalar, numWellDofs, numWellDofs>>;
    using WVector = Dune::BlockVector<Dune::FieldVector<Scalar, numWellDofs>>;

    virtual void addBCDMatrix(std::vector<BMatrix>& b_matrices,
        std::vector<CMatrix>& c_matrices,
        std::vector<DMatrix>& d_matrices,
        Opm::SparseTable<int>& wcells) const = 0;

protected:
    // simulation parameters
    std::vector<RateVector> connectionRates_;
    std::vector<Scalar> B_avg_;
    bool changed_to_stopped_this_step_ = false;
    bool thp_update_iterations = false;
    int number_of_well_reopenings_{0};

    Scalar wpolymer() const;
    Scalar wfoam() const;
    Scalar wsalt() const;
    Scalar wmicrobes() const;
    Scalar woxygen() const;
    Scalar wurea() const;

    virtual Scalar getRefDensity() const = 0;

    std::vector<Scalar>
    initialWellRateFractions(const Simulator& ebosSimulator,
                             const WellStateType& well_state) const;

    // check whether the well is operable under BHP limit with current reservoir condition
    virtual void checkOperabilityUnderBHPLimit(const WellStateType& well_state,
                                               const Simulator& simulator,
                                               DeferredLogger& deferred_logger) = 0;

    // check whether the well is operable under THP limit with current reservoir condition
    virtual void checkOperabilityUnderTHPLimit(const Simulator& simulator,
                                               const WellStateType& well_state,
                                               const GroupStateHelperType& groupStateHelper) = 0;

    virtual void updateIPR(const Simulator& simulator,
                           DeferredLogger& deferred_logger) const = 0;

    virtual void assembleWellEqWithoutIteration(const Simulator& simulator,
                                                const GroupStateHelperType& groupStateHelper,
                                                const double dt,
                                                const WellInjectionControls& inj_controls,
                                                const WellProductionControls& prod_controls,
                                                WellStateType& well_state,
                                                const bool solving_with_zero_rate) = 0;

    // iterate well equations with the specified control until converged
    virtual bool iterateWellEqWithControl(const Simulator& simulator,
                                          const double dt,
                                          const WellInjectionControls& inj_controls,
                                          const WellProductionControls& prod_controls,
                                          const GroupStateHelperType& groupStateHelper,
                                          WellStateType& well_state) = 0;

    bool iterateWellEquations(const Simulator& simulator,
                              const double dt,
                              const GroupStateHelperType& groupStateHelper,
                              WellStateType& well_state);

    bool solveWellWithOperabilityCheck(const Simulator& simulator,
                                       const double dt,
                                       const Well::InjectionControls& inj_controls,
                                       const Well::ProductionControls& prod_controls,
                                       const GroupStateHelperType& groupStateHelper,
                                       WellStateType& well_state);

    std::optional<Scalar>
    estimateOperableBhp(const Simulator& ebos_simulator,
                        const double dt,
                        const GroupStateHelperType& groupStateHelper,
                        const SummaryState& summary_state,
                        WellStateType& well_state);

    bool solveWellWithBhp(const Simulator& simulator,
                          const double dt,
                          const Scalar bhp,
                          const GroupStateHelperType& groupStateHelper,
                          WellStateType& well_state);

    bool solveWellWithZeroRate(const Simulator& simulator,
                               const double dt,
                               const GroupStateHelperType& groupStateHelper,
                               WellStateType& well_state);

    bool solveWellForTesting(const Simulator& simulator,
                             const GroupStateHelperType& groupStateHelper,
                             WellStateType& well_state);


    template<class GasLiftSingleWell>
    std::unique_ptr<GasLiftSingleWell> initializeGliftWellTest_(const Simulator& simulator,
                                                                WellStateType& well_state,
                                                                const GroupState<Scalar>& group_state,
                                                                GLiftEclWells& ecl_well_map,
                                                                DeferredLogger& deferred_logger);

    Eval getPerfCellPressure(const FluidState& fs) const;

    // get the transmissibility multiplier for specific perforation
    template<class Value, class Callback>
    void getTransMult(Value& trans_mult,
                      const Simulator& simulator,
                      const int cell_idx,
                      Callback& extendEval) const;

    // get the well transmissibility for specific perforation
    template<class Value>
    void getTw(std::vector<Value>&         wi,
               const int                   perf,
               const IntensiveQuantities&  intQuants,
               const Value&                trans_mult,
               const SingleWellStateType&  ws) const;

    // get the mobility for specific perforation
    template<class Value, class Callback>
    void getMobility(const Simulator& simulator,
                     const int local_perf_index,
                     std::vector<Value>& mob,
                     Callback& extendEval,
                     [[maybe_unused]] DeferredLogger& deferred_logger) const;

    void computeConnLevelProdInd(const FluidState& fs,
                                 const std::function<Scalar(const Scalar)>& connPICalc,
                                 const std::vector<Scalar>& mobility,
                                 Scalar* connPI) const;

    void computeConnLevelInjInd(const FluidState& fs,
                                const Phase preferred_phase,
                                const std::function<Scalar(const Scalar)>& connIICalc,
                                const std::vector<Scalar>& mobility,
                                Scalar* connII,
                                DeferredLogger& deferred_logger) const;

    Scalar computeConnectionDFactor(const int perf,
                                    const IntensiveQuantities& intQuants,
                                    const SingleWellStateType& ws) const;
};

} // namespace Opm

#include "WellInterface_impl.hpp"

#endif // OPM_WELLINTERFACE_HEADER_INCLUDED