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;
                            const Scalar oilPV = pv_cell * fs.saturation(FluidSystem::oilPhaseIdx).value();
                            const Scalar gasPV = pv_cell * fs.saturation(FluidSystem::gasPhaseIdx).value();
                            attr.rs_oil += fs.Rs().value() * oilPV;
                            attr.oil_pv += oilPV;
                            attr.rv_gas += fs.Rv().value() * gasPV;
                            attr.gas_pv += gasPV;
                        }
                        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])
                    , rs_oil(data[8])
                    , oil_pv(data[9])
                    , rv_gas(data[10])
                    , gas_pv(data[11])
                {}
 
                Attributes(const Attributes& rhs)
                    : Attributes()
                {
                    this->data = rhs.data;
                }
 
                Attributes& operator=(const Attributes& rhs)
                {
                    this->data = rhs.data;
                    return *this;
                }
 
                std::array<Scalar,12> data;
                Scalar& pressure;
                Scalar& temperature;
                Scalar& rs;
                Scalar& rv;
                Scalar& rsw;
                Scalar& rvw;
                Scalar& pv;
                Scalar& saltConcentration;
                // Oil (gas) volume weighted Rs (Rv) and corresponding weights.
                Scalar& rs_oil;
                Scalar& oil_pv;
                Scalar& rv_gas;
                Scalar& gas_pv;
            };
 
            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 */