diff --git a/ebos/ecltransmissibility.hh b/ebos/ecltransmissibility.hh index c880e6422..9c714b2d5 100644 --- a/ebos/ecltransmissibility.hh +++ b/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 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 keywords {"TRANX", "TRANY", "TRANZ"}; + std::array,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& is_tran, + const std::array,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(); const auto& elemEndIt = gridView.template end(); @@ -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()); - - if (c1 > c2) + 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); - 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]; - else - // Scale transmissibilities with scale factor from inputTranx - trans_[isId] *= inputTranxData[c1]; + trans_[isId] = trans[0][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]; - else - // Scale transmissibilities with scale factor from inputTrany - trans_[isId] *= inputTranyData[c1]; + trans_[isId] = trans[1][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]; - else - // Scale transmissibilities with scale factor from inputTranz - trans_[isId] *= inputTranzData[c1]; + trans_[isId] = trans[2][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,3> + createTransmissibilityArrays(const std::array& 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,3> trans = + { std::vector(is_tran[0] ? numElem : 0, 0), + std::vector(is_tran[1] ? numElem : 0, 0), + std::vector(is_tran[2] ? numElem : 0, 0)}; + + // compute the transmissibilities for all intersections + auto elemIt = gridView.template begin(); + const auto& elemEndIt = gridView.template end(); + + 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 void computeFaceProperties(const Intersection& intersection, diff --git a/opm/simulators/utils/ParallelEclipseState.cpp b/opm/simulators/utils/ParallelEclipseState.cpp index d11a57235..477b65012 100644 --- a/opm/simulators/utils/ParallelEclipseState.cpp +++ b/opm/simulators/utils/ParallelEclipseState.cpp @@ -68,17 +68,18 @@ const std::vector& 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>&>(m_intProps)[keyword]; - local_data.resize(m_activeSize()); + auto& local_data = const_cast>&>(m_intProps)[keyword]; + 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 ParallelFieldPropsManager::get_global_int(const std::string& keyword) const @@ -121,16 +122,17 @@ const std::vector& 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>&>(m_doubleProps)[keyword]; - local_data.resize(m_activeSize()); + auto& local_data = const_cast>&>(m_doubleProps)[keyword]; + 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 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& 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& buffer) +{ + Opm::deserialize_tran(m_tran, buffer); +} bool ParallelFieldPropsManager::has_double(const std::string& keyword) const { diff --git a/opm/simulators/utils/ParallelEclipseState.hpp b/opm/simulators/utils/ParallelEclipseState.hpp index 60a189ce1..e0f4ffe40 100644 --- a/opm/simulators/utils/ParallelEclipseState.hpp +++ b/opm/simulators/utils/ParallelEclipseState.hpp @@ -20,6 +20,7 @@ #define PARALLEL_ECLIPSE_STATE_HPP #include +#include #include #include @@ -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& trans) const override; + + void deserialize_tran(const std::vector& buffer) override; protected: - std::map> m_intProps; //!< Map of integer properties in process-local compressed indices. - std::map> m_doubleProps; //!< Map of double properties in process-local compressed indices. + std::map> m_intProps; //!< Map of integer properties in process-local compressed indices. + std::map> 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 m_comm; //!< Collective communication handler. std::function m_activeSize; //!< active size function of the grid std::function m_local2Global; //!< mapping from local to global cartesian indices + std::unordered_map m_tran; //!< calculators map }; diff --git a/opm/simulators/utils/PropsCentroidsDataHandle.hpp b/opm/simulators/utils/PropsCentroidsDataHandle.hpp index fbbea7bb7..14875fd41 100644 --- a/opm/simulators/utils/PropsCentroidsDataHandle.hpp +++ b/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; //! \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 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(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(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(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(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(data->second[counter++]); + { + const auto& pair = data->second[counter++]; + m_distributed_fieldProps.m_intProps[intKey].data[index] = static_cast(pair.first); + m_distributed_fieldProps.m_intProps[intKey].value_status[index] = static_cast(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(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 @@ -216,7 +253,9 @@ private: //! \brief The names of the keys of the double fields. std::vector m_doubleKeys; /// \brief The data per element as a vector mapped from the local id. - std::unordered_map > elementData_; + /// + /// each entry is a pair of data and value_status. + std::unordered_map > > elementData_; /// \brief The cell centroids of the distributed grid. std::vector& m_centroids; /// \brief The amount of data to send for each element