RateConverter.hpp

Go to the documentation of this file.

/*
  Copyright 2014, 2015 SINTEF ICT, Applied Mathematics.
  Copyright 2014, 2015 Statoil ASA.
  Copyright 2017, 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/>.
*/

#ifndef OPM_RATECONVERTER_HPP_HEADER_INCLUDED
#define OPM_RATECONVERTER_HPP_HEADER_INCLUDED

#include <opm/grid/utility/RegionMapping.hpp>

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

#include <opm/simulators/utils/DeferredLoggingErrorHelpers.hpp>
#include <opm/simulators/utils/ParallelCommunication.hpp>

#include <dune/grid/common/gridenums.hh>
#include <dune/grid/common/rangegenerators.hh>

#include <array>
#include <cassert>
#include <unordered_map>
#include <utility>
#include <vector>

namespace Opm {
    namespace RateConverter {

        template <class FluidSystem, class Region>
        class SurfaceToReservoirVoidage {
        public:
            using Scalar = typename FluidSystem::Scalar;
            SurfaceToReservoirVoidage(const Region&     region)
                : rmap_      (region)
                , attr_      (rmap_, Attributes())
            {}

            template <typename ElementContext, class Simulator>
            void defineState(const Simulator& simulator)
            {
                // create map from cell to region and set all attributes to
                // zero
                for (const auto& reg : rmap_.activeRegions()) {
                    auto& ra = attr_.attributes(reg);
                    ra.pressure = 0.0;
                    ra.temperature = 0.0;
                    ra.rs = 0.0;
                    ra.rv = 0.0;
                    ra.pv = 0.0;
                    ra.saltConcentration = 0.0;
                    ra.rsw = 0.0;
                    ra.rvw = 0.0;
                }

                // quantities for pore volume average
                std::unordered_map<RegionId, Attributes> attributes_pv;

                // quantities for hydrocarbon volume average
                std::unordered_map<RegionId, Attributes> attributes_hpv;

                for (const auto& reg : rmap_.activeRegions()) {
                    attributes_pv.insert({reg, Attributes()});
                    attributes_hpv.insert({reg, Attributes()});
                }

                ElementContext elemCtx( simulator );
                const auto& gridView = simulator.gridView();
                const auto& comm = gridView.comm();

                OPM_BEGIN_PARALLEL_TRY_CATCH();
                for (const auto& elem : elements(gridView, Dune::Partitions::interior)) {
                    elemCtx.updatePrimaryStencil(elem);
                    elemCtx.updatePrimaryIntensiveQuantities(/*timeIdx=*/0);
                    const unsigned cellIdx = elemCtx.globalSpaceIndex(/*spaceIdx=*/0, /*timeIdx=*/0);
                    const auto& intQuants = elemCtx.intensiveQuantities(/*spaceIdx=*/0, /*timeIdx=*/0);
                    const auto& fs = intQuants.fluidState();
                    // use pore volume weighted averages.
                    const Scalar pv_cell =
                            simulator.model().dofTotalVolume(simulator.vanguard().gridEquilIdxToGridIdx(cellIdx))
                            * intQuants.porosity().value();

                    // only count oil and gas filled parts of the domain
                    Scalar hydrocarbon = 1.0;
                    if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
                        hydrocarbon -= fs.saturation(FluidSystem::waterPhaseIdx).value();
                    }

                    const int reg = rmap_.region(cellIdx);
                    assert(reg >= 0);

                    // sum p, rs, rv, and T.
                    const Scalar hydrocarbonPV = pv_cell*hydrocarbon;
                    if (hydrocarbonPV > 0.) {
                        auto& attr = attributes_hpv[reg];
                        attr.pv += hydrocarbonPV;
                        if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
                            attr.rs += fs.Rs().value() * hydrocarbonPV;
                            attr.rv += fs.Rv().value() * hydrocarbonPV;
                        }
                        if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
                            attr.pressure += fs.pressure(FluidSystem::oilPhaseIdx).value() * hydrocarbonPV;
                            attr.temperature += fs.temperature(FluidSystem::oilPhaseIdx).value() * hydrocarbonPV;
                        } else {
                            assert(FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx));
                            attr.pressure += fs.pressure(FluidSystem::gasPhaseIdx).value() * hydrocarbonPV;
                            attr.temperature += fs.temperature(FluidSystem::gasPhaseIdx).value() * hydrocarbonPV;
                        }
                        attr.saltConcentration += fs.saltConcentration().value() * hydrocarbonPV;
                        if (FluidSystem::enableDissolvedGasInWater()) {
                            attr.rsw += fs.Rsw().value() * hydrocarbonPV; // scale with total volume?
                        }
                        if (FluidSystem::enableVaporizedWater()) {
                            attr.rvw += fs.Rvw().value() * hydrocarbonPV; // scale with total volume?
                        }
                    }

                    if (pv_cell > 0.) {
                        auto& attr = attributes_pv[reg];
                        attr.pv += pv_cell;
                        if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
                            attr.rs += fs.Rs().value() * pv_cell;
                            attr.rv += fs.Rv().value() * pv_cell;
                        }
                        if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
                            attr.pressure += fs.pressure(FluidSystem::oilPhaseIdx).value() * pv_cell;
                            attr.temperature += fs.temperature(FluidSystem::oilPhaseIdx).value() * pv_cell;
                        } else if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
                             attr.pressure += fs.pressure(FluidSystem::gasPhaseIdx).value() * pv_cell;
                             attr.temperature += fs.temperature(FluidSystem::gasPhaseIdx).value() * pv_cell;
                        } else {
                            assert(FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx));
                            attr.pressure += fs.pressure(FluidSystem::waterPhaseIdx).value() * pv_cell;
                            attr.temperature += fs.temperature(FluidSystem::waterPhaseIdx).value() * pv_cell;
                        }
                        attr.saltConcentration += fs.saltConcentration().value() * pv_cell;
                        if (FluidSystem::enableDissolvedGasInWater()) {
                            attr.rsw += fs.Rsw().value() * pv_cell;
                        }
                        if (FluidSystem::enableVaporizedWater()) {
                            attr.rvw += fs.Rvw().value() * pv_cell;
                        }
                    }
                }

                OPM_END_PARALLEL_TRY_CATCH("SurfaceToReservoirVoidage::defineState() failed: ", simulator.vanguard().grid().comm());

                this->sumRates(attributes_hpv,
                               attributes_pv,
                               comm);
            }

            using RegionId = typename RegionMapping<Region>::RegionId;

            template <class Coeff>
            void
            calcCoeff(const RegionId r, const int pvtRegionIdx, Coeff& coeff) const;

            template <class Coeff , class Rates>
            void
            calcCoeff(const RegionId r, const int pvtRegionIdx, const Rates& surface_rates, Coeff& coeff) const;

            template <class Coeff>
            void
            calcCoeff(  const int     pvtRegionIdx,
                        const Scalar  p,
                        const Scalar  rs,
                        const Scalar  rv,
                        const Scalar  rsw,
                        const Scalar  rvw,
                        const Scalar  T,
                        const Scalar  saltConcentration,
                        Coeff&        coeff) const;

            template <class Coeff>
            void
            calcInjCoeff(const RegionId r, const int pvtRegionIdx, Coeff& coeff) const;

            template <class Rates>
            void calcReservoirVoidageRates(const RegionId r,
                                           const int      pvtRegionIdx,
                                           const Rates&   surface_rates,
                                           Rates&         voidage_rates) const;

            template <typename SurfaceRates, typename VoidageRates>
            void calcReservoirVoidageRates(const int           pvtRegionIdx,
                                           const Scalar        p,
                                           const Scalar        rs,
                                           const Scalar        rv,
                                           const Scalar        rsw,
                                           const Scalar        rvw,
                                           const Scalar        T,
                                           const Scalar        saltConcentration,
                                           const SurfaceRates& surface_rates,
                                           VoidageRates&       voidage_rates) const;

            template <class Rates>
            std::pair<Scalar, Scalar>
            inferDissolvedVaporisedRatio(const Scalar rsMax,
                                         const Scalar rvMax,
                                         const Rates& surface_rates) const;

            template <class SolventModule>
            void
            calcCoeffSolvent(const RegionId r, const int pvtRegionIdx, Scalar& coeff) const
            {
                const auto& ra = attr_.attributes(r);
                const Scalar p = ra.pressure;
                const Scalar T = ra.temperature;
                const auto& solventPvt = SolventModule::solventPvt();
                const Scalar bs = solventPvt.inverseFormationVolumeFactor(pvtRegionIdx, T, p);
                coeff = 1.0 / bs;
            }

        private:
            const RegionMapping<Region> rmap_;

            struct Attributes {
                Attributes()
                    : data{0.0}
                    , pressure(data[0])
                    , temperature(data[1])
                    , rs(data[2])
                    , rv(data[3])
                    , rsw(data[4])
                    , rvw(data[5])
                    , pv(data[6])
                    , saltConcentration(data[7])
                {}

                Attributes(const Attributes& rhs)
                    : Attributes()
                {
                    this->data = rhs.data;
                }

                Attributes& operator=(const Attributes& rhs)
                {
                    this->data = rhs.data;
                    return *this;
                }

                std::array<Scalar,8> data;
                Scalar& pressure;
                Scalar& temperature;
                Scalar& rs;
                Scalar& rv;
                Scalar& rsw;
                Scalar& rvw;
                Scalar& pv;
                Scalar& saltConcentration;
            };

            void sumRates(std::unordered_map<RegionId,Attributes>& attributes_hpv,
                          std::unordered_map<RegionId,Attributes>& attributes_pv,
                          Parallel::Communication comm);

            RegionAttributeHelpers::RegionAttributes<RegionId, Attributes> attr_;
        };

    } // namespace RateConverter
} // namespace Opm

#endif  /* OPM_RATECONVERTER_HPP_HEADER_INCLUDED */