opm-simulators/opm/simulators/wells/WellGroupConstraints.cpp

/*
  Copyright 2017 SINTEF Digital, Mathematics and Cybernetics.
  Copyright 2017 Statoil ASA.
  Copyright 2018 IRIS

  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/>.
*/

#include <config.h>

#include <opm/simulators/wells/WellGroupConstraints.hpp>

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

#include <opm/material/fluidsystems/BlackOilDefaultFluidSystemIndices.hpp>

#include <opm/simulators/wells/GroupStateHelper.hpp>
#include <opm/simulators/wells/WellInterfaceGeneric.hpp>
#include <opm/simulators/wells/WellState.hpp>

namespace Opm
{

template<typename Scalar, typename IndexTraits>
std::pair<bool, Scalar>
WellGroupConstraints<Scalar, IndexTraits>::
checkGroupConstraintsInj(const Group& group,
                         const GroupStateHelperType& groupStateHelper,
                         const Scalar efficiencyFactor,
                         const SummaryState& summaryState,
                         const RateConvFunc& rateConverter,
                         const bool check_guide_rate) const
{
    const auto& well_state = groupStateHelper.wellState();

    // Translate injector type from control to Phase.
    const auto& well_controls = well_.wellEcl().injectionControls(summaryState);
    auto injectorType = well_controls.injector_type;
    Phase injectionPhase;
    switch (injectorType) {
    case InjectorType::WATER:
    {
        injectionPhase = Phase::WATER;
        break;
    }
    case InjectorType::OIL:
    {
        injectionPhase = Phase::OIL;
        break;
    }
    case InjectorType::GAS:
    {
        injectionPhase = Phase::GAS;
        break;
    }
    default:
        throw("Expected WATER, OIL or GAS as type for injector " + well_.name());
    }

    // Make conversion factors for RESV <-> surface rates.
    std::vector<Scalar> resv_coeff(well_.phaseUsage().numActivePhases(), 1.0);
    rateConverter(0, well_.pvtRegionIdx(), group.name(), resv_coeff); // FIPNUM region 0 here, should use FIPNUM from WELSPECS.

    const auto& ws = well_state.well(well_.indexOfWell());
    // Call check for the well's injection phase.
    return groupStateHelper.checkGroupConstraintsInj(
        well_.name(),
        well_.wellEcl().groupName(),
        group,
        ws.surface_rates.data(),
        injectionPhase,
        efficiencyFactor,
        resv_coeff,
        check_guide_rate
    );
}

template<typename Scalar, typename IndexTraits>
std::pair<bool, Scalar>
WellGroupConstraints<Scalar, IndexTraits>::
checkGroupConstraintsProd(const Group& group,
                          const GroupStateHelperType& groupStateHelper,
                          const Scalar efficiencyFactor,
                          const RateConvFunc& rateConverter,
                          const bool check_guide_rate) const
{
    const auto& well_state = groupStateHelper.wellState();

    // Make conversion factors for RESV <-> surface rates.
    std::vector<Scalar> resv_coeff(well_.phaseUsage().numActivePhases(), 1.0);
    rateConverter(0, well_.pvtRegionIdx(), group.name(), resv_coeff); // FIPNUM region 0 here, should use FIPNUM from WELSPECS.

    const auto& ws = well_state.well(well_.indexOfWell());
    return groupStateHelper.checkGroupConstraintsProd(well_.name(),
        well_.wellEcl().groupName(),
        group,
        ws.surface_rates.data(),
        efficiencyFactor,
        resv_coeff,
        check_guide_rate
    );
}

template<typename Scalar, typename IndexTraits>
bool WellGroupConstraints<Scalar, IndexTraits>::
checkGroupConstraints(const GroupStateHelperType& groupStateHelper,
                      const Schedule& schedule,
                      const SummaryState& summaryState,
                      const RateConvFunc& rateConverter,
                      const bool check_guide_rate,
                      WellStateType& well_state) const
{
    const auto& well = well_.wellEcl();
    const int well_index = well_.indexOfWell();
    auto& ws = well_state.well(well_index);

    if (well.isInjector()) {
        const auto currentControl = ws.injection_cmode;

        if (currentControl != Well::InjectorCMode::GRUP) {
            // This checks only the first encountered group limit,
            // in theory there could be several, and then we should
            // test all but the one currently applied. At that point,
            // this if-statement should be removed and we should always
            // check, skipping over only the single group parent whose
            // control is the active one for the well (if any).
            const auto& group = schedule.getGroup(well.groupName(), well_.currentStep());
            const Scalar efficiencyFactor = well.getEfficiencyFactor() *
                                            well_state[well.name()].efficiency_scaling_factor;
            const std::pair<bool, Scalar> group_constraint =
                this->checkGroupConstraintsInj(group, groupStateHelper, efficiencyFactor,
                                               summaryState,
                                               rateConverter,
                                               check_guide_rate);
            // If a group constraint was broken, we set the current well control to
            // be GRUP.
            if (group_constraint.first) {
                ws.injection_cmode = Well::InjectorCMode::GRUP;
                const int np = well_state.numPhases();
                for (int p = 0; p<np; ++p) {
                    ws.surface_rates[p] *= group_constraint.second;
                }
            }
            return group_constraint.first;
        }
    }

    if (well.isProducer( )) {
        const auto currentControl = ws.production_cmode;

        if (currentControl != Well::ProducerCMode::GRUP) {
            // This checks only the first encountered group limit,
            // in theory there could be several, and then we should
            // test all but the one currently applied. At that point,
            // this if-statement should be removed and we should always
            // check, skipping over only the single group parent whose
            // control is the active one for the well (if any).
            const auto& group = schedule.getGroup(well.groupName(), well_.currentStep());
            const Scalar efficiencyFactor = well.getEfficiencyFactor() *
                                            well_state[well.name()].efficiency_scaling_factor;
            const std::pair<bool, Scalar> group_constraint =
                this->checkGroupConstraintsProd(group, groupStateHelper,
                                                efficiencyFactor,
                                                rateConverter,
                                                check_guide_rate);
            // If a group constraint was broken, we set the current well control to
            // be GRUP.
            if (group_constraint.first) {
                ws.production_cmode = Well::ProducerCMode::GRUP;
                const int np = well_state.numPhases();
                for (int p = 0; p<np; ++p) {
                    ws.surface_rates[p] *= group_constraint.second;
                }
            }
            return group_constraint.first;
        }
    }

    return false;
}

template class WellGroupConstraints<double, BlackOilDefaultFluidSystemIndices>;

#if FLOW_INSTANTIATE_FLOAT
template class WellGroupConstraints<float, BlackOilDefaultFluidSystemIndices>;
#endif

} // namespace Opm