WellBhpThpCalculator.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_WELL_BPH_THP_CALCULATOR_HEADER_INCLUDED
#define OPM_WELL_BPH_THP_CALCULATOR_HEADER_INCLUDED

#include <functional>
#include <optional>
#include <vector>

namespace Opm {

class DeferredLogger;
class SummaryState;
class Well;
template<typename Scalar, typename IndexTraits> class WellInterfaceGeneric;
template<typename Scalar, typename IndexTraits> class WellState;

template<typename Scalar, typename IndexTraits>
class WellBhpThpCalculator {
public:
    explicit WellBhpThpCalculator(const WellInterfaceGeneric<Scalar, IndexTraits>& well) : well_(well) {}

    bool wellHasTHPConstraints(const SummaryState& summaryState) const;

    Scalar getTHPConstraint(const SummaryState& summaryState) const;

    Scalar mostStrictBhpFromBhpLimits(const SummaryState& summaryState) const;

    Scalar calculateThpFromBhp(const std::vector<Scalar>& rates,
                               const Scalar bhp,
                               const Scalar rho,
                               const std::optional<Scalar>& alq,
                               const Scalar thp_limit,
                               DeferredLogger& deferred_logger) const;

    std::optional<Scalar>
    computeBhpAtThpLimitProd(const std::function<std::vector<Scalar>(const Scalar)>& frates,
                             const SummaryState& summary_state,
                             const Scalar maxPerfPress,
                             const Scalar rho,
                             const Scalar alq_value,
                             const Scalar thp_limit,
                             DeferredLogger& deferred_logger) const;

    std::optional<Scalar>
    computeBhpAtThpLimitInj(const std::function<std::vector<Scalar>(const Scalar)>& frates,
                            const SummaryState& summary_state,
                            const Scalar rho,
                            const Scalar flo_rel_tol,
                            const int max_iteration,
                            const bool throwOnError,
                            DeferredLogger& deferred_logger) const;

    void updateThp(const Scalar rho,
                   const std::function<Scalar()>& alq_value,
                   WellState<Scalar, IndexTraits>& well_state,
                   const SummaryState& summary_state,
                   DeferredLogger& deferred_logger) const;

    template<class EvalWell>
    EvalWell calculateBhpFromThp(const WellState<Scalar, IndexTraits>& well_state,
                                 const std::vector<EvalWell>& rates,
                                 const Well& well,
                                 const SummaryState& summaryState,
                                 const Scalar rho,
                                 DeferredLogger& deferred_logger) const;

    Scalar calculateMinimumBhpFromThp(const WellState<Scalar, IndexTraits>& well_state,
                                      const Well& well,
                                      const SummaryState& summaryState,
                                      const Scalar rho) const;

  bool isStableSolution(const WellState<Scalar, IndexTraits>& well_state,
                        const Well& well,
                        const std::vector<Scalar>& rates,
                        const SummaryState& summaryState) const;

  std::optional<Scalar>
  estimateStableBhp (const WellState<Scalar, IndexTraits>& well_state,
                    const Well& well,
                    const std::vector<Scalar>& rates,
                    const Scalar rho,
                    const SummaryState& summaryState) const;

  std::optional<Scalar>
  estimateStableBhp (const WellState<Scalar, IndexTraits>& well_state,
                    const Well& well,
                    const std::vector<Scalar>& rates,
                    const Scalar rho,
                    const SummaryState& summaryState,
                    const Scalar alq) const;

  std::pair<Scalar, Scalar>
  getFloIPR(const WellState<Scalar, IndexTraits>& well_state,
            const Well& well,
            const SummaryState& summary_state) const;

  static bool bruteForceBracketCommonTHP(const std::function<Scalar(const Scalar)>& eq,
                                const std::array<Scalar, 2>& range,
                                Scalar& low, Scalar& high,
                                std::optional<Scalar>& approximate_solution,
                                const Scalar& limit,
                                DeferredLogger& deferred_logger);

  static bool bruteForceBracketCommonTHP(const std::function<Scalar(const Scalar)>& eq,
                                Scalar& min_thp, Scalar& max_thp);

private:
    template<class ErrorPolicy>
    std::optional<Scalar>
    computeBhpAtThpLimitInjImpl(const std::function<std::vector<Scalar>(const Scalar)>& frates,
                                const SummaryState& summary_state,
                                const Scalar rho,
                                const Scalar flo_rel_tol,
                                const int max_iteration,
                                DeferredLogger& deferred_logger) const;

    std::optional<Scalar>
    bhpMax(const std::function<Scalar(const Scalar)>& fflo,
           const Scalar bhp_limit,
           const Scalar maxPerfPress,
           const Scalar vfp_flo_front,
           DeferredLogger& deferred_logger) const;

    std::optional<Scalar>
    computeBhpAtThpLimit(const std::function<std::vector<Scalar>(const Scalar)>& frates,
                         const std::function<Scalar(const std::vector<Scalar>)>& fbhp,
                         const std::array<Scalar, 2>& range,
                         DeferredLogger& deferred_logger) const;

    Scalar getVfpBhpAdjustment(const Scalar bph_tab, const Scalar thp_limit) const;

    bool bisectBracket(const std::function<Scalar(const Scalar)>& eq,
                       const std::array<Scalar, 2>& range,
                       Scalar& low, Scalar& high,
                       std::optional<Scalar>& approximate_solution,
                       DeferredLogger& deferred_logger) const;

  static bool bruteForceBracket(const std::function<Scalar(const Scalar)>& eq,
                                const std::array<Scalar, 2>& range,
                                Scalar& low, Scalar& high,
                                DeferredLogger& deferred_logger);


    Scalar findThpFromBhpIteratively(const std::function<Scalar(const Scalar, const Scalar)>& thp_func,
                                     const Scalar bhp,
                                     const Scalar thp_limit,
                                     const Scalar dp,
                                     DeferredLogger& deferred_logger) const;

    const WellInterfaceGeneric<Scalar, IndexTraits>& well_; 
};

}

#endif // OPM_WELL_BHP_THP_CALCULATOR_HEADER_INCLUDED