Merge pull request #2793 from blattms/fix-tran-application

Fixes application of transmissibility modifiers from the edit section.

ebos/ecltransmissibility.hh

@@ -389,7 +389,7 @@ public:
         }
 
         // potentially overwrite and/or modify  transmissibilities based on input from deck
-        updateFromEclState_();
+        updateFromEclState_(global);
 
         // Create mapping from global to local index
         const size_t cartesianSize = cartMapper.cartesianSize();
@@ -512,20 +512,48 @@ private:
         }
     }
 
-    void updateFromEclState_()
+    void updateFromEclState_(bool global)
+    {
+        const FieldPropsManager* fp =
+            (global) ? &(vanguard_.eclState().fieldProps()) :
+            &(vanguard_.eclState().globalFieldProps());
+
+        std::array<bool,3> is_tran {fp->tran_active("TRANX"),
+                                    fp->tran_active("TRANY"),
+                                    fp->tran_active("TRANZ")};
+
+        if( !(is_tran[0] ||is_tran[1] || is_tran[2]) )
+        {
+            // Skip unneeded expensive traversals
+            return;
+        }
+
+        std::array<std::string, 3> keywords {"TRANX", "TRANY", "TRANZ"};
+        std::array<std::vector<double>,3> trans = createTransmissibilityArrays(is_tran);
+        auto key = keywords.begin();
+        auto perform = is_tran.begin();
+
+        for (auto it = trans.begin(); it != trans.end(); ++it, ++key, ++perform)
+        {
+            if(perform)
+                fp->apply_tran(*key, *it);
+        }
+
+        resetTransmissibilityFromArrays(is_tran, trans);
+    }
+
+    /// \brief overwrites calculated transmissibilities
+    ///
+    /// \param is_tran Whether TRAN{XYZ} have been modified.
+    /// \param trans Arrays with modified transmissibilities TRAN{XYZ}
+    void resetTransmissibilityFromArrays(const std::array<bool,3>& is_tran,
+                                         const std::array<std::vector<double>,3>& trans)
     {
         const auto& gridView = vanguard_.gridView();
         const auto& cartMapper = vanguard_.cartesianIndexMapper();
         const auto& cartDims = cartMapper.cartesianDimensions();
         ElementMapper elemMapper(gridView, Dune::mcmgElementLayout());
 
-        const auto& fp = vanguard_.eclState().fieldProps();
-        const auto& inputTranxData = fp.get_double("TRANX");
-        const auto& inputTranyData = fp.get_double("TRANY");
-        const auto& inputTranzData = fp.get_double("TRANZ");
-        bool tranx_deckAssigned = false;                     // Ohh my ....
-        bool trany_deckAssigned = false;
-        bool tranz_deckAssigned = false;
         // compute the transmissibilities for all intersections
         auto elemIt = gridView.template begin</*codim=*/ 0>();
         const auto& elemEndIt = gridView.template end</*codim=*/ 0>();
@@ -540,46 +568,113 @@ private:
                 if (!intersection.neighbor())
                     continue; // intersection is on the domain boundary
 
+                // In the EclState TRANX[c1] is transmissibility in X+
+                // direction. Ordering of compressed (c1,c2) and cartesian index
+                // (gc1, gc2) is coherent (c1 < c2 <=> gc1 < gc2). This also
+                // holds for the global grid. While distributing changes the
+                // order of the local indices, the transmissibilities are still
+                // stored at the cell with the lower global cartesian index as
+                // the fieldprops are communicated by the grid.
                 unsigned c1 = elemMapper.index(intersection.inside());
                 unsigned c2 = elemMapper.index(intersection.outside());
-
+                int gc1 = cartMapper.cartesianIndex(c1);
-                if (c1 > c2)
+                int gc2 = cartMapper.cartesianIndex(c2);
+                if (gc1 > gc2)
                     continue; // we only need to handle each connection once, thank you.
 
                 auto isId = isId_(c1, c2);
 
-                int gc1 = std::min(cartMapper.cartesianIndex(c1), cartMapper.cartesianIndex(c2));
-                int gc2 = std::max(cartMapper.cartesianIndex(c1), cartMapper.cartesianIndex(c2));
-
                 if (gc2 - gc1 == 1 && cartDims[0] > 1) {
-                    if (tranx_deckAssigned)
+                    if (is_tran[0])
                         // set simulator internal transmissibilities to values from inputTranx
-                        trans_[isId] = inputTranxData[c1];
+                        trans_[isId] = trans[0][c1];
-                    else
-                        // Scale transmissibilities with scale factor from inputTranx
-                        trans_[isId] *= inputTranxData[c1];
                 }
                 else if (gc2 - gc1 == cartDims[0] && cartDims[1] > 1) {
-                    if (trany_deckAssigned)
+                    if (is_tran[1])
                         // set simulator internal transmissibilities to values from inputTrany
-                        trans_[isId] = inputTranyData[c1];
+                        trans_[isId] = trans[1][c1];
-                    else
-                        // Scale transmissibilities with scale factor from inputTrany
-                        trans_[isId] *= inputTranyData[c1];
                 }
                 else if (gc2 - gc1 == cartDims[0]*cartDims[1]) {
-                    if (tranz_deckAssigned)
+                    if (is_tran[2])
                         // set simulator internal transmissibilities to values from inputTranz
-                        trans_[isId] = inputTranzData[c1];
+                        trans_[isId] = trans[2][c1];
-                    else
-                        // Scale transmissibilities with scale factor from inputTranz
-                        trans_[isId] *= inputTranzData[c1];
                 }
                 //else.. We don't support modification of NNC at the moment.
             }
         }
     }
 
+    /// \brief Creates TRANS{XYZ} arrays for modification by FieldProps data
+    ///
+    /// \param is_tran Whether TRAN{XYZ} will be modified. If entry is false the array will be empty
+    /// \returns an array of vector (TRANX, TRANY, TRANZ}
+    std::array<std::vector<double>,3>
+    createTransmissibilityArrays(const std::array<bool,3>& is_tran)
+    {
+        const auto& gridView = vanguard_.gridView();
+        const auto& cartMapper = vanguard_.cartesianIndexMapper();
+        const auto& cartDims = cartMapper.cartesianDimensions();
+        ElementMapper elemMapper(gridView, Dune::mcmgElementLayout());
+
+        auto numElem = vanguard_.gridView().size(/*codim=*/0);
+        std::array<std::vector<double>,3> trans =
+            { std::vector<double>(is_tran[0] ? numElem : 0, 0),
+              std::vector<double>(is_tran[1] ? numElem : 0, 0),
+              std::vector<double>(is_tran[2] ? numElem : 0, 0)};
+
+        // compute the transmissibilities for all intersections
+        auto elemIt = gridView.template begin</*codim=*/ 0>();
+        const auto& elemEndIt = gridView.template end</*codim=*/ 0>();
+
+        for (; elemIt != elemEndIt; ++elemIt) {
+            const auto& elem = *elemIt;
+            auto isIt = gridView.ibegin(elem);
+            const auto& isEndIt = gridView.iend(elem);
+            for (; isIt != isEndIt; ++ isIt) {
+                // store intersection, this might be costly
+                const auto& intersection = *isIt;
+                if (!intersection.neighbor())
+                    continue; // intersection is on the domain boundary
+
+                // In the EclState TRANX[c1] is transmissibility in X+
+                // direction. Ordering of compressed (c1,c2) and cartesian index
+                // (gc1, gc2) is coherent (c1 < c2 <=> gc1 < gc2). This also
+                // holds for the global grid. While distributing changes the
+                // order of the local indices, the transmissibilities are still
+                // stored at the cell with the lower global cartesian index as
+                // the fieldprops are communicated by the grid.
+                unsigned c1 = elemMapper.index(intersection.inside());
+                unsigned c2 = elemMapper.index(intersection.outside());
+                int gc1 = cartMapper.cartesianIndex(c1);
+                int gc2 = cartMapper.cartesianIndex(c2);
+
+                if (gc1 > gc2)
+                    continue; // we only need to handle each connection once, thank you.
+
+                auto isId = isId_(c1, c2);
+
+                if (gc2 - gc1 == 1 && cartDims[0] > 1) {
+                    if (is_tran[0])
+                        // set simulator internal transmissibilities to values from inputTranx
+                         trans[0][c1] = trans_[isId];
+                }
+                else if (gc2 - gc1 == cartDims[0] && cartDims[1] > 1) {
+                    if (is_tran[1])
+                        // set simulator internal transmissibilities to values from inputTrany
+                         trans[1][c1] = trans_[isId];
+                }
+                else if (gc2 - gc1 == cartDims[0]*cartDims[1]) {
+                    if (is_tran[2])
+                        // set simulator internal transmissibilities to values from inputTranz
+                         trans[2][c1] = trans_[isId];
+                }
+                //else.. We don't support modification of NNC at the moment.
+            }
+        }
+
+        return trans;
+    }
+
 
     template <class Intersection>
     void computeFaceProperties(const Intersection& intersection,

opm/simulators/utils/ParallelEclipseState.cpp

@@ -68,17 +68,18 @@ const std::vector<int>& ParallelFieldPropsManager::get_int(const std::string& ke
         // Some of the keywords might be defaulted.
         // We will let rank 0 create them and distribute them using get_global_int
         auto data = get_global_int(keyword);
-        auto& local_data = const_cast<std::map<std::string, std::vector<int>>&>(m_intProps)[keyword];
+        auto& local_data = const_cast<std::map<std::string, Fieldprops::FieldData<int>>&>(m_intProps)[keyword];
-        local_data.resize(m_activeSize());
+        local_data.data.resize(m_activeSize());
+        local_data.value_status.resize(m_activeSize());
 
         for (int i = 0; i < m_activeSize(); ++i)
         {
-            local_data[i] = data[m_local2Global(i)];
+            local_data.data[i] = data[m_local2Global(i)];
         }
-        return local_data;
+        return local_data.data;
     }
 
-    return it->second;
+    return it->second.data;
 }
 
 std::vector<int> ParallelFieldPropsManager::get_global_int(const std::string& keyword) const
@@ -121,16 +122,17 @@ const std::vector<double>& ParallelFieldPropsManager::get_double(const std::stri
         // Some of the keywords might be defaulted.
         // We will let rank 0 create them and distribute them using get_global_int
         auto data = get_global_double(keyword);
-        auto& local_data = const_cast<std::map<std::string, std::vector<double>>&>(m_doubleProps)[keyword];
+        auto& local_data = const_cast<std::map<std::string, Fieldprops::FieldData<double>>&>(m_doubleProps)[keyword];
-        local_data.resize(m_activeSize());
+        local_data.data.resize(m_activeSize());
+        local_data.value_status.resize(m_activeSize());
         for (int i = 0; i < m_activeSize(); ++i)
         {
-            local_data[i] = data[m_local2Global(i)];
+            local_data.data[i] = data[m_local2Global(i)];
         }
-        return local_data;
+        return local_data.data;
     }
 
-    return it->second;
+    return it->second.data;
 }
 
 
@@ -165,6 +167,17 @@ std::vector<double> ParallelFieldPropsManager::get_global_double(const std::stri
     return result;
 }
 
+bool ParallelFieldPropsManager::tran_active(const std::string& keyword) const
+{
+    auto calculator = m_tran.find(keyword);
+    return calculator != m_tran.end() && calculator->second.size();
+}
+
+void ParallelFieldPropsManager::apply_tran(const std::string& keyword,
+                                           std::vector<double>& data) const
+{
+    Opm::apply_tran(m_tran, m_doubleProps, m_activeSize(), keyword, data);
+}
 
 bool ParallelFieldPropsManager::has_int(const std::string& keyword) const
 {
@@ -172,6 +185,10 @@ bool ParallelFieldPropsManager::has_int(const std::string& keyword) const
     return it != m_intProps.end();
 }
 
+void ParallelFieldPropsManager::deserialize_tran(const std::vector<char>& buffer)
+{
+    Opm::deserialize_tran(m_tran, buffer);
+}
 
 bool ParallelFieldPropsManager::has_double(const std::string& keyword) const
 {

opm/simulators/utils/ParallelEclipseState.hpp

@@ -20,6 +20,7 @@
 #define PARALLEL_ECLIPSE_STATE_HPP
 
 #include <opm/parser/eclipse/EclipseState/EclipseState.hpp>
+#include <opm/parser/eclipse/EclipseState/Grid/TranCalculator.hpp>
 #include <dune/common/parallel/mpihelper.hh>
 
 #include <functional>
@@ -95,13 +96,20 @@ public:
         m_local2Global = std::bind(&T::cartesianIndex, mapper,
                                    std::placeholders::_1);
     }
+
+    bool tran_active(const std::string& keyword) const override;
+
+    void apply_tran(const std::string& keyword, std::vector<double>& trans) const override;
+
+    void deserialize_tran(const std::vector<char>& buffer) override;
 protected:
-    std::map<std::string, std::vector<int>> m_intProps; //!< Map of integer properties in process-local compressed indices.
+    std::map<std::string, Opm::Fieldprops::FieldData<int>> m_intProps; //!< Map of integer properties in process-local compressed indices.
-    std::map<std::string, std::vector<double>> m_doubleProps; //!< Map of double properties in process-local compressed indices.
+    std::map<std::string, Opm::Fieldprops::FieldData<double>> m_doubleProps; //!< Map of double properties in process-local compressed indices.
     FieldPropsManager& m_manager; //!< Underlying field property manager (only used on root process).
     Dune::CollectiveCommunication<Dune::MPIHelper::MPICommunicator> m_comm; //!< Collective communication handler.
     std::function<int(void)> m_activeSize; //!< active size function of the grid
     std::function<int(const int)> m_local2Global; //!< mapping from local to global cartesian indices
+    std::unordered_map<std::string, Fieldprops::TranCalculator> m_tran; //!< calculators map
 };
 
 

opm/simulators/utils/PropsCentroidsDataHandle.hpp

@@ -50,7 +50,7 @@ class PropsCentroidsDataHandle
 {
 public:
     //! \brief the data type we send (ints are converted to double)
-    using DataType = double;
+    using DataType = std::pair<double, unsigned char>;
 
     //! \brief Constructor
     //! \param grid The grid where the loadbalancing is happening.
@@ -81,9 +81,18 @@ public:
             int position = 0;
             Mpi::pack(m_intKeys, buffer, position, comm);
             Mpi::pack(m_doubleKeys, buffer, position, comm);
+            int calcStart = position;
+            {
+                std::vector<char> tran_buffer = globalProps.serialize_tran();
+                position += tran_buffer.size();
+                buffer.insert(buffer.end(), std::make_move_iterator(tran_buffer.begin()), std::make_move_iterator(tran_buffer.end()));
+            }
             comm.broadcast(&position, 1, 0);
             comm.broadcast(buffer.data(), position, 0);
 
+            // Unpack Calculator as we need it here, too.
+            m_distributed_fieldProps.deserialize_tran( std::vector<char>(buffer.begin() + calcStart, buffer.end()) );
+
             // copy data to persistent map based on local id
             m_no_data = m_intKeys.size() + m_doubleKeys.size() +
                 Grid::dimensionworld;
@@ -101,15 +110,26 @@ public:
                 data.reserve(m_no_data);
 
                 for (const auto& intKey : m_intKeys)
-                    data.push_back(globalProps.get_int(intKey)[index]);
+                {
+                    const auto& fieldData = globalProps.get_int_field_data(intKey);
+                    data.emplace_back(fieldData.data[index],
+                                      static_cast<unsigned char>(fieldData.value_status[index]));
+                }
 
                 for (const auto& doubleKey : m_doubleKeys)
-                    data.push_back(globalProps.get_double(doubleKey)[index]);
+                {
+                    // We need to allow unsupported keywords to get the data
+                    // for TranCalculator, too.
+                    const auto& fieldData = globalProps.get_double_field_data(doubleKey,
+                                                                              /* allow_unsupported = */ true);
+                    data.emplace_back(fieldData.data[index],
+                                      static_cast<unsigned char>(fieldData.value_status[index]));
+                }
 
                 auto cartIndex = cartMapper.cartesianIndex(index);
                 const auto& center = eclGridOnRoot->getCellCenter(cartIndex);
                 for (int dim = 0; dim < Grid::dimensionworld; ++dim)
-                    data.push_back(center[dim]);
+                    data.emplace_back(center[dim], '1'); // write garbage for value_status
             }
         }
         else
@@ -121,6 +141,7 @@ public:
             int position{};
             Mpi::unpack(m_intKeys, buffer, position, comm);
             Mpi::unpack(m_doubleKeys, buffer, position, comm);
+            m_distributed_fieldProps.deserialize_tran( std::vector<char>(buffer.begin() + position, buffer.end()) );
             m_no_data = m_intKeys.size() + m_doubleKeys.size() +
                 Grid::dimensionworld;
         }
@@ -130,10 +151,16 @@ public:
     {
         // distributed grid is now correctly set up.
         for(const auto& intKey : m_intKeys)
-            m_distributed_fieldProps.m_intProps[intKey].resize(m_grid.size(0));
+        {
+            m_distributed_fieldProps.m_intProps[intKey].data.resize(m_grid.size(0));
+            m_distributed_fieldProps.m_intProps[intKey].value_status.resize(m_grid.size(0));
+        }
 
         for(const auto& doubleKey : m_doubleKeys)
-            m_distributed_fieldProps.m_doubleProps[doubleKey].resize(m_grid.size(0));
+        {
+            m_distributed_fieldProps.m_doubleProps[doubleKey].data.resize(m_grid.size(0));
+            m_distributed_fieldProps.m_doubleProps[doubleKey].value_status.resize(m_grid.size(0));
+        }
 
         m_centroids.resize(m_grid.size(0) * Grid::dimensionworld);
 
@@ -153,15 +180,23 @@ public:
             assert(data != elementData_.end());
 
             for(const auto& intKey : m_intKeys)
-                m_distributed_fieldProps.m_intProps[intKey][index] = static_cast<int>(data->second[counter++]);
+            {
+                const auto& pair = data->second[counter++];
+                m_distributed_fieldProps.m_intProps[intKey].data[index] = static_cast<int>(pair.first);
+                m_distributed_fieldProps.m_intProps[intKey].value_status[index] = static_cast<value::status>(pair.second);
+            }
 
             for(const auto& doubleKey : m_doubleKeys)
-                m_distributed_fieldProps.m_doubleProps[doubleKey][index] = data->second[counter++];
+            {
+                const auto& pair = data->second[counter++];
+                m_distributed_fieldProps.m_doubleProps[doubleKey].data[index] = pair.first;
+                m_distributed_fieldProps.m_doubleProps[doubleKey].value_status[index] = static_cast<value::status>(pair.second);
+            }
 
             auto centroidIter = m_centroids.begin() + Grid::dimensionworld * index;
             auto centroidIterEnd = centroidIter + Grid::dimensionworld;
             for ( ; centroidIter != centroidIterEnd; ++centroidIter )
-                *centroidIter = data->second[counter++];
+                *centroidIter = data->second[counter++].first; // value_status discarded
         }
     }
 
@@ -191,7 +226,9 @@ public:
         auto iter = elementData_.find(m_grid.localIdSet().id(e));
         assert(iter != elementData_.end());
         for(const auto& data : iter->second)
+        {
             buffer.write(data);
+        }
     }
 
     template<class BufferType, class EntityType>
@@ -216,7 +253,9 @@ private:
     //! \brief The names of the keys of the double fields.
     std::vector<std::string> m_doubleKeys;
     /// \brief The data per element as a vector mapped from the local id.
-    std::unordered_map<typename LocalIdSet::IdType, std::vector<double> > elementData_;
+    ///
+    /// each entry is a pair of data and value_status.
+    std::unordered_map<typename LocalIdSet::IdType, std::vector<std::pair<double,unsigned char> > > elementData_;
     /// \brief The cell centroids of the distributed grid.
     std::vector<double>& m_centroids;
     /// \brief The amount of data to send for each element