RegionAttributeHelpers.hpp

Go to the documentation of this file.

/*
  Copyright 2014, 2015 SINTEF ICT, Applied Mathematics.
  Copyright 2014, 2015 Statoil ASA.
  Copyright 2017, IRIS
  Copyright 2017, Equinor
 
  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_REGIONATTRIBUTEHELPERS_HPP_HEADER_INCLUDED
#define OPM_REGIONATTRIBUTEHELPERS_HPP_HEADER_INCLUDED
 
#include <opm/simulators/utils/BlackoilPhases.hpp>
 
#include <opm/grid/utility/RegionMapping.hpp>
 
#include <dune/grid/common/gridenums.hh>
 
#include <algorithm>
#include <memory>
#include <stdexcept>
#include <type_traits>
#include <unordered_map>
#include <utility>
#include <vector>
 
namespace Opm {
    namespace RegionAttributeHelpers {
            namespace Select {
                template <class RegionID, bool>
                struct RegionIDParameter
                {
                    using type =
                        typename std::remove_reference<RegionID>::type &;
                };
 
                template <class RegionID>
                struct RegionIDParameter<RegionID, true>
                {
                    using type = RegionID;
                };
            } // Select
 
            template<class Scalar, bool is_parallel>
            struct AverageIncrementCalculator
            {
                std::tuple<Scalar, Scalar, Scalar, Scalar, int>
                operator()(const std::vector<Scalar>& pressure,
                           const std::vector<Scalar>& temperature,
                           const std::vector<Scalar>& rs,
                           const std::vector<Scalar>& rv,
                           const std::vector<Scalar>& ownership,
                           std::size_t cell){
                    if ( ownership[cell] )
                    {
                        return std::make_tuple(pressure[cell],
                                               temperature[cell],
                                               rs[cell],
                                               rv[cell],
                                               1);
                    }
                    else
                    {
                        return std::make_tuple(0, 0, 0, 0, 0);
                    }
                }
            };
            template<class Scalar>
            struct AverageIncrementCalculator<Scalar, false>
            {
                std::tuple<Scalar, Scalar, Scalar, Scalar, int>
                operator()(const std::vector<Scalar>& pressure,
                           const std::vector<Scalar>& temperature,
                           const std::vector<Scalar>& rs,
                           const std::vector<Scalar>& rv,
                           const std::vector<Scalar>&,
                           std::size_t cell){
                    return std::make_tuple(pressure[cell],
                                           temperature[cell],
                                           rs[cell],
                                           rv[cell],
                                           1);
                }
            };
            template <typename RegionId, class Attributes>
            class RegionAttributes
            {
            public:
                using RegionID =
                    typename Select::RegionIDParameter
                    <RegionId, std::is_integral<RegionId>::value>::type;
 
                using ID =
                    typename std::remove_reference<RegionId>::type;
 
                struct Value {
                    Value(const Attributes& attr)
                        : attr_(attr)
                        , cell_(-1)
                    {}
 
                    Attributes attr_;
                    int        cell_;
                };
 
                using AttributeMap =
                    std::unordered_map<ID, std::unique_ptr<Value>>;
 
 
                template <class RMap>
                RegionAttributes(const RMap&       rmap,
                                 const Attributes& attr)
                {
                    using VT = typename AttributeMap::value_type;
 
                    for (const auto& r : rmap.activeRegions()) {
                        auto v = std::make_unique<Value>(attr);
 
                        const auto stat = attr_.insert(VT(r, std::move(v)));
 
                        if (stat.second) {
                            // New value inserted.
                            const auto& cells = rmap.cells(r);
 
                            assert (! cells.empty());
 
                            // Region's representative cell.
                            stat.first->second->cell_ = cells[0];
                        }
                    }
                }
 
                int cell(const RegionID reg) const
                {
                    return this->find(reg).cell_;
                }
 
                bool has(const RegionID reg) const
                {
                    return this->attr_.find(reg) != this->attr_.end();
                }
 
                void insert(const RegionID r, const Attributes& attr)
                {
                    auto [pos, inserted] = this->attr_.try_emplace(r, std::make_unique<Value>(attr));
                    if (inserted) {
                        pos->second->cell_ = -1; // NOT -1.0 -- "cell_" is 'int'
                    }
                }
 
                const AttributeMap& attributes() const
                {
                    return attr_;
                }
 
 
                const Attributes& attributes(const RegionID reg) const
                {
                    return this->find(reg).attr_;
                }
 
                Attributes& attributes(const RegionID reg)
                {
                    return this->find(reg).attr_;
                }
 
            private:
 
                AttributeMap attr_;
 
                const Value& find(const RegionID reg) const
                {
                    const auto& i = attr_.find(reg);
 
                    if (i == attr_.end()) {
                        throw std::invalid_argument("Unknown region ID");
                    }
 
                    return *i->second;
                }
 
                Value& find(const RegionID reg)
                {
                    const auto& i = attr_.find(reg);
 
                    if (i == attr_.end()) {
                        throw std::invalid_argument("Unknown region ID");
                    }
 
                    return *i->second;
                }
            };
 
            namespace PhaseUsed {
                inline bool
                water(const PhaseUsage& pu)
                {
                    return pu.phase_used[ BlackoilPhases::Aqua ] != 0;
                }
 
                inline bool
                oil(const PhaseUsage& pu)
                {
                    return pu.phase_used[ BlackoilPhases::Liquid ] != 0;
                }
 
                inline bool
                gas(const PhaseUsage& pu)
                {
                    return pu.phase_used[ BlackoilPhases::Vapour ] != 0;
                }
            } // namespace PhaseUsed
 
            namespace PhasePos {
                inline int
                water(const PhaseUsage& pu)
                {
                    int p = -1;
 
                    if (PhaseUsed::water(pu)) {
                        p = pu.phase_pos[ BlackoilPhases::Aqua ];
                    }
 
                    return p;
                }
 
                inline int
                oil(const PhaseUsage& pu)
                {
                    int p = -1;
 
                    if (PhaseUsed::oil(pu)) {
                        p = pu.phase_pos[ BlackoilPhases::Liquid ];
                    }
 
                    return p;
                }
 
                inline int
                gas(const PhaseUsage& pu)
                {
                    int p = -1;
 
                    if (PhaseUsed::gas(pu)) {
                        p = pu.phase_pos[ BlackoilPhases::Vapour ];
                    }
 
                    return p;
                }
            } // namespace PhasePos
 
    } // namespace RegionAttributesHelpers
} // namespace Opm
 
#endif  /* OPM_REGIONATTRIBUTEHELPERS_HPP_HEADER_INCLUDED */