diff --git a/CMakeLists.txt b/CMakeLists.txt index ff8bb9c34..4750f7645 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -317,8 +317,8 @@ endmacro (files_hook) macro (tests_hook) endmacro (tests_hook) - - + + # all setup common to the OPM library modules is done here include (OpmLibMain) @@ -559,12 +559,6 @@ foreach(OBJ ${FLOW_VARIANT_MODELS}) set_property(TARGET flow_lib${OBJ} PROPERTY EXCLUDE_FROM_ALL ${FLOW_VARIANTS_DEFAULT_ENABLE_IF}) endforeach() -add_library(flow_libblackoil_poly OBJECT flow/flow_ebos_blackoil_poly.cpp) -if(TARGET fmt::fmt) - target_link_libraries(flow_libblackoil_poly fmt::fmt) - endif() -set_property(TARGET flow_libblackoil_poly PROPERTY EXCLUDE_FROM_ALL ${FLOW_POLY_ONLY_DEFAULT_ENABLE_IF}) - opm_add_test(flow ONLY_COMPILE ALWAYS_ENABLE @@ -577,17 +571,15 @@ opm_add_test(flow $ ) -opm_add_test(flow_poly +opm_add_test(flow_blackoil_polyhedralgrid ONLY_COMPILE ALWAYS_ENABLE DEFAULT_ENABLE_IF ${FLOW_POLY_ONLY_DEFAULT_ENABLE_IF} DEPENDS opmsimulators LIBRARIES opmsimulators SOURCES - flow/flow_blackoil_poly.cpp - $ + flow/flow_blackoil_polyhedralgrid.cpp $) -target_compile_definitions(flow_poly PRIVATE USE_POLYHEDRALGRID) opm_add_test(flow_distribute_z ONLY_COMPILE @@ -608,7 +600,7 @@ if(dune-alugrid_FOUND) set(FLOW_ALUGRID_ONLY_DEFAULT_ENABLE_IF "TRUE") endif() - opm_add_test(flow_alugrid + opm_add_test(flow_blackoil_alugrid ONLY_COMPILE ALWAYS_ENABLE DEFAULT_ENABLE_IF ${FLOW_ALUGRID_ONLY_DEFAULT_ENABLE_IF} @@ -617,7 +609,6 @@ if(dune-alugrid_FOUND) SOURCES flow/flow_blackoil_alugrid.cpp $) - target_compile_definitions(flow_alugrid PRIVATE USE_ALUGRID) endif() if (BUILD_FLOW) @@ -685,5 +676,3 @@ endif() install(DIRECTORY doc/man1 DESTINATION ${CMAKE_INSTALL_MANDIR} FILES_MATCHING PATTERN "*.1") - - diff --git a/ebos/collecttoiorank.cc b/ebos/collecttoiorank.cc index 08b39af9b..dc60d529e 100644 --- a/ebos/collecttoiorank.cc +++ b/ebos/collecttoiorank.cc @@ -23,1121 +23,18 @@ #include #include +#include -#include #include -#include -#include -#include -#include #if HAVE_DUNE_FEM #include #include #include #endif // HAVE_DUNE_FEM -#if HAVE_DUNE_ALUGRID -#include -#include -#include "alucartesianindexmapper.hh" -#endif // HAVE_DUNE_ALUGRID - - -#include -#include -#include -#include -#include - -namespace { - std::vector toVector(const std::set& string_set) - { - return { string_set.begin(), string_set.end() }; - } -} - namespace Opm { -// global id -class GlobalCellIndex -{ - int globalId_; - int localIndex_; - bool isInterior_; - -public: - GlobalCellIndex() - : globalId_(-1) - , localIndex_(-1) - , isInterior_(true) - {} - void setGhost() - { isInterior_ = false; } - - void setId(int globalId) - { globalId_ = globalId; } - void setIndex(int localIndex) - { localIndex_ = localIndex; } - - int localIndex () const - { return localIndex_; } - int id () const - { return globalId_; } - bool isInterior() const - { return isInterior_; } -}; - -using IndexMapType = std::vector; -using IndexMapStorageType = std::vector; -using P2PCommunicatorType = Dune::Point2PointCommunicator; -using MessageBufferType = typename P2PCommunicatorType::MessageBufferType; -class DistributeIndexMapping : public P2PCommunicatorType::DataHandleInterface -{ -protected: - const std::vector& distributedGlobalIndex_; - IndexMapType& localIndexMap_; - IndexMapStorageType& indexMaps_; - std::map globalPosition_; - std::set& recv_; - std::vector& ranks_; - -public: - DistributeIndexMapping(const std::vector& globalIndex, - const std::vector& distributedGlobalIndex, - IndexMapType& localIndexMap, - IndexMapStorageType& indexMaps, - std::vector& ranks, - std::set& recv, - bool isIORank) - : distributedGlobalIndex_(distributedGlobalIndex) - , localIndexMap_(localIndexMap) - , indexMaps_(indexMaps) - , globalPosition_() - , recv_(recv) - , ranks_(ranks) - { - size_t size = globalIndex.size(); - // create mapping globalIndex --> localIndex - if ( isIORank ) // ioRank - for (size_t index = 0; index < size; ++index) - globalPosition_.insert(std::make_pair(globalIndex[index], index)); - - // we need to create a mapping from local to global - if (!indexMaps_.empty()) { - if (isIORank) - ranks_.resize(size, -1); - // for the ioRank create a localIndex to index in global state map - IndexMapType& indexMap = indexMaps_.back(); - size_t localSize = localIndexMap_.size(); - indexMap.resize(localSize); - for (size_t i=0; isecond; - // Using max should be backwards compatible - ranks_[entry] = std::max(ranks_[entry], rank); - } - if (rankIt != recv_.end()) - ++rankIt; - } - #ifndef NDEBUG - for (const auto& rank: ranks_) - assert(rank>=0); - #endif - } - } - - void pack(int link, MessageBufferType& buffer) - { - // we should only get one link - if (link != 0) - throw std::logic_error("link in method pack is not 0 as execpted"); - - // pack all interior global cell id's - int size = localIndexMap_.size(); - buffer.write(size); - - for (int index = 0; index < size; ++index) { - int globalIdx = distributedGlobalIndex_[localIndexMap_[index]]; - buffer.write(globalIdx); - } - } - - void unpack(int link, MessageBufferType& buffer) - { - // get index map for current link - IndexMapType& indexMap = indexMaps_[link]; - assert(!globalPosition_.empty()); - - // unpack all interior global cell id's - int numCells = 0; - buffer.read(numCells); - indexMap.resize(numCells); - for (int index = 0; index < numCells; ++index) { - buffer.read(indexMap[index]); - } - } -}; - - -/// \brief Communication handle to scatter the global index -template -class ElementIndexScatterHandle -{ -public: - ElementIndexScatterHandle(const EquilMapper& sendMapper, const Mapper& recvMapper, std::vector& elementIndices) - : sendMapper_(sendMapper), recvMapper_(recvMapper), elementIndices_(elementIndices) - {} - using DataType = int; - bool fixedSize(int /*dim*/, int /*codim*/) - { - return true; - } - - template - std::size_t size(const T&) - { - return 1; - } - template - void gather(B& buffer, const T& t) - { - buffer.write(sendMapper_.index(t)); - } - template - void scatter(B& buffer, const T& t, std::size_t) - { - buffer.read(elementIndices_[recvMapper_.index(t)]); - } - - bool contains(int dim, int codim) - { - return dim==3 && codim==0; - } -private: - const EquilMapper& sendMapper_; - const Mapper& recvMapper_; - std::vector& elementIndices_; -}; - -/// \brief Communication handle to scatter the global index -template -class ElementIndexHandle -{ -public: - ElementIndexHandle(const Mapper& mapper, std::vector& elementIndices) - : mapper_(mapper), elementIndices_(elementIndices) - {} - using DataType = int; - bool fixedSize(int /*dim*/, int /*codim*/) - { - return true; - } - - template - std::size_t size(const T&) - { - return 1; - } - template - void gather(B& buffer, const T& t) - { - buffer.write(elementIndices_[mapper_.index(t)]); - } - template - void scatter(B& buffer, const T& t, std::size_t) - { - buffer.read(elementIndices_[mapper_.index(t)]); - } - - bool contains(int dim, int codim) - { - return dim==3 && codim==0; - } -private: - const Mapper& mapper_; - std::vector& elementIndices_; -}; - -class PackUnPackCellData : public P2PCommunicatorType::DataHandleInterface -{ - const data::Solution& localCellData_; - data::Solution& globalCellData_; - - const IndexMapType& localIndexMap_; - const IndexMapStorageType& indexMaps_; - -public: - PackUnPackCellData(const data::Solution& localCellData, - data::Solution& globalCellData, - const IndexMapType& localIndexMap, - const IndexMapStorageType& indexMaps, - size_t globalSize, - bool isIORank) - : localCellData_(localCellData) - , globalCellData_(globalCellData) - , localIndexMap_(localIndexMap) - , indexMaps_(indexMaps) - { - if (isIORank) { - // add missing data to global cell data - for (const auto& pair : localCellData_) { - const std::string& key = pair.first; - std::size_t containerSize = globalSize; - [[maybe_unused]] auto ret = globalCellData_.insert(key, pair.second.dim, - std::vector(containerSize), - pair.second.target); - assert(ret.second); - } - - MessageBufferType buffer; - pack(0, buffer); - - // the last index map is the local one - doUnpack(indexMaps.back(), buffer); - } - } - - // pack all data associated with link - void pack(int link, MessageBufferType& buffer) - { - // we should only get one link - if (link != 0) - throw std::logic_error("link in method pack is not 0 as expected"); - - // write all cell data registered in local state - for (const auto& pair : localCellData_) { - const auto& data = pair.second.data; - - // write all data from local data to buffer - write(buffer, localIndexMap_, data); - } - } - - void doUnpack(const IndexMapType& indexMap, MessageBufferType& buffer) - { - // we loop over the data as - // its order governs the order the data got received. - for (auto& pair : localCellData_) { - const std::string& key = pair.first; - auto& data = globalCellData_.data(key); - - //write all data from local cell data to buffer - read(buffer, indexMap, data); - } - } - - // unpack all data associated with link - void unpack(int link, MessageBufferType& buffer) - { doUnpack(indexMaps_[link], buffer); } - -protected: - template - void write(MessageBufferType& buffer, - const IndexMapType& localIndexMap, - const Vector& vector, - unsigned int offset = 0, - unsigned int stride = 1) const - { - unsigned int size = localIndexMap.size(); - buffer.write(size); - assert(vector.size() >= stride * size); - for (unsigned int i=0; i - void read(MessageBufferType& buffer, - const IndexMapType& indexMap, - Vector& vector, - unsigned int offset = 0, - unsigned int stride = 1) const - { - unsigned int size = 0; - buffer.read(size); - assert(size == indexMap.size()); - for (unsigned int i=0; ipack(0, buffer); - - // pass a dummy_link to satisfy virtual class - const int dummyLink = -1; - this->unpack(dummyLink, buffer); - } - - // pack all data associated with link - void pack(int link, MessageBufferType& buffer) - { - // we should only get one link - if (link != 0) { - throw std::logic_error { - "link in method pack is not 0 as expected" - }; - } - - // write all group and network (node/branch) data - this->localGroupAndNetworkData_.write(buffer); - } - - // unpack all data associated with link - void unpack(int /*link*/, MessageBufferType& buffer) - { this->globalGroupAndNetworkData_.read(buffer); } -}; - -class PackUnPackBlockData : public P2PCommunicatorType::DataHandleInterface -{ - const std::map, double>& localBlockData_; - std::map, double>& globalBlockValues_; - -public: - PackUnPackBlockData(const std::map, double>& localBlockData, - std::map, double>& globalBlockValues, - bool isIORank) - : localBlockData_(localBlockData) - , globalBlockValues_(globalBlockValues) - { - if (isIORank) { - MessageBufferType buffer; - pack(0, buffer); - - // pass a dummyLink to satisfy virtual class - int dummyLink = -1; - unpack(dummyLink, buffer); - } - } - - // pack all data associated with link - void pack(int link, MessageBufferType& buffer) - { - // we should only get one link - if (link != 0) - throw std::logic_error("link in method pack is not 0 as expected"); - - // write all block data - unsigned int size = localBlockData_.size(); - buffer.write(size); - for (const auto& map : localBlockData_) { - buffer.write(map.first.first); - buffer.write(map.first.second); - buffer.write(map.second); - } - } - - // unpack all data associated with link - void unpack(int /*link*/, MessageBufferType& buffer) - { - // read all block data - unsigned int size = 0; - buffer.read(size); - for (size_t i = 0; i < size; ++i) { - std::string name; - int idx; - double data; - buffer.read(name); - buffer.read(idx); - buffer.read(data); - globalBlockValues_[std::make_pair(name, idx)] = data; - } - } -}; - -class PackUnPackWBPData : public P2PCommunicatorType::DataHandleInterface -{ - const data::WellBlockAveragePressures& localWBPData_; - data::WellBlockAveragePressures& globalWBPValues_; - -public: - PackUnPackWBPData(const data::WellBlockAveragePressures& localWBPData, - data::WellBlockAveragePressures& globalWBPValues, - const bool isIORank) - : localWBPData_ (localWBPData) - , globalWBPValues_(globalWBPValues) - { - if (! isIORank) { - return; - } - - MessageBufferType buffer; - pack(0, buffer); - - // Pass a dummy link to satisfy base class API requirement - const int dummyLink = -1; - unpack(dummyLink, buffer); - } - - // Pack all data associated with link - void pack(int link, MessageBufferType& buffer) - { - // We should only get one link - if (link != 0) { - throw std::logic_error { - "link (" + std::to_string(link) + - ") in method pack() is not 0 as expected" - }; - } - - // Write all local, per-well, WBP data - this->localWBPData_.write(buffer); - } - - // Unpack all data associated with link - void unpack([[maybe_unused]] const int link, - MessageBufferType& buffer) - { - this->globalWBPValues_.read(buffer); - } -}; - -class PackUnPackWellTestState : public P2PCommunicatorType::DataHandleInterface -{ -public: - PackUnPackWellTestState(const WellTestState& localWTestState, - WellTestState& globalWTestState, - bool isIORank) - : local_(localWTestState) - , global_(globalWTestState) - { - if (isIORank) { - MessageBufferType buffer; - pack(0, buffer); - - // pass a dummyLink to satisfy virtual class - int dummyLink = -1; - unpack(dummyLink, buffer); - } - } - - void pack(int link, MessageBufferType& buffer) { - if (link != 0) - throw std::logic_error("link in method pack is not 0 as expected"); - this->local_.pack(buffer); - } - - void unpack(int, MessageBufferType& buffer) { - this->global_.unpack(buffer); - } - -private: - const WellTestState& local_; - WellTestState& global_; -}; - -class PackUnPackAquiferData : public P2PCommunicatorType::DataHandleInterface -{ - const data::Aquifers& localAquiferData_; - data::Aquifers& globalAquiferData_; - -public: - PackUnPackAquiferData(const data::Aquifers& localAquiferData, - data::Aquifers& globalAquiferData, - bool isIORank) - : localAquiferData_(localAquiferData) - , globalAquiferData_(globalAquiferData) - { - if (isIORank) { - MessageBufferType buffer; - pack(0, buffer); - - // pass a dummyLink to satisfy virtual class - int dummyLink = -1; - unpack(dummyLink, buffer); - } - } - - // pack all data associated with link - void pack(int link, MessageBufferType& buffer) - { - // we should only get one link - if (link != 0) - throw std::logic_error("link in method pack is not 0 as expected"); - - int size = localAquiferData_.size(); - buffer.write(size); - for (const auto& [key, data] : localAquiferData_) { - buffer.write(key); - data.write(buffer); - } - } - - // unpack all data associated with link - void unpack(int /*link*/, MessageBufferType& buffer) - { - int size; - buffer.read(size); - for (int i = 0; i < size; ++i) { - int key; - buffer.read(key); - data::AquiferData data; - data.read(buffer); - - auto& aq = this->globalAquiferData_[key]; - - this->unpackCommon(data, aq); - - if (auto const* aquFet = data.typeData.get(); - aquFet != nullptr) - { - this->unpackAquFet(*aquFet, aq); - } - else if (auto const* aquCT = data.typeData.get(); - aquCT != nullptr) - { - this->unpackCarterTracy(*aquCT, aq); - } - else if (auto const* aquNum = data.typeData.get(); - aquNum != nullptr) - { - this->unpackNumericAquifer(*aquNum, aq); - } - } - } - -private: - void unpackCommon(const data::AquiferData& data, data::AquiferData& aq) - { - aq.aquiferID = std::max(aq.aquiferID, data.aquiferID); - aq.pressure = std::max(aq.pressure, data.pressure); - aq.initPressure = std::max(aq.initPressure, data.initPressure); - aq.datumDepth = std::max(aq.datumDepth, data.datumDepth); - aq.fluxRate += data.fluxRate; - aq.volume += data.volume; - } - - void unpackAquFet(const data::FetkovichData& aquFet, data::AquiferData& aq) - { - if (! aq.typeData.is()) { - auto* tData = aq.typeData.create(); - *tData = aquFet; - } - else { - const auto& src = aquFet; - auto& dst = *aq.typeData.getMutable(); - - dst.initVolume = std::max(dst.initVolume , src.initVolume); - dst.prodIndex = std::max(dst.prodIndex , src.prodIndex); - dst.timeConstant = std::max(dst.timeConstant, src.timeConstant); - } - } - - void unpackCarterTracy(const data::CarterTracyData& aquCT, data::AquiferData& aq) - { - if (! aq.typeData.is()) { - auto* tData = aq.typeData.create(); - *tData = aquCT; - } - else { - const auto& src = aquCT; - auto& dst = *aq.typeData.getMutable(); - - dst.timeConstant = std::max(dst.timeConstant , src.timeConstant); - dst.influxConstant = std::max(dst.influxConstant, src.influxConstant); - dst.waterDensity = std::max(dst.waterDensity , src.waterDensity); - dst.waterViscosity = std::max(dst.waterViscosity, src.waterViscosity); - - dst.dimensionless_time = std::max(dst.dimensionless_time , src.dimensionless_time); - dst.dimensionless_pressure = std::max(dst.dimensionless_pressure, src.dimensionless_pressure); - } - } - - void unpackNumericAquifer(const data::NumericAquiferData& aquNum, data::AquiferData& aq) - { - if (! aq.typeData.is()) { - auto* tData = aq.typeData.create(); - *tData = aquNum; - } - else { - const auto& src = aquNum; - auto& dst = *aq.typeData.getMutable(); - - std::transform(src.initPressure.begin(), - src.initPressure.end(), - dst.initPressure.begin(), - dst.initPressure.begin(), - [](const double p0_1, const double p0_2) - { - return std::max(p0_1, p0_2); - }); - } - } -}; - -class PackUnpackInterRegFlows : public P2PCommunicatorType::DataHandleInterface -{ - const EclInterRegFlowMap& localInterRegFlows_; - EclInterRegFlowMap& globalInterRegFlows_; - -public: - PackUnpackInterRegFlows(const EclInterRegFlowMap& localInterRegFlows, - EclInterRegFlowMap& globalInterRegFlows, - const bool isIORank) - : localInterRegFlows_(localInterRegFlows) - , globalInterRegFlows_(globalInterRegFlows) - { - if (! isIORank) { return; } - - MessageBufferType buffer; - this->pack(0, buffer); - - // pass a dummy_link to satisfy virtual class - const int dummyLink = -1; - this->unpack(dummyLink, buffer); - } - - // pack all data associated with link - void pack(int link, MessageBufferType& buffer) - { - // we should only get one link - if (link != 0) { - throw std::logic_error { - "link in method pack is not 0 as expected" - }; - } - - // write all inter-region flow data - this->localInterRegFlows_.write(buffer); - } - - // unpack all data associated with link - void unpack(int /*link*/, MessageBufferType& buffer) - { this->globalInterRegFlows_.read(buffer); } -}; - -class PackUnpackFlows : public P2PCommunicatorType::DataHandleInterface -{ - const std::array, std::vector>>, 3>& localFlows_; - std::array, std::vector>>, 3>& globalFlows_; - -public: - PackUnpackFlows(const std::array, std::vector>>, 3> & localFlows, - std::array, std::vector>>, 3>& globalFlows, - const bool isIORank) - : localFlows_(localFlows) - , globalFlows_(globalFlows) - { - if (! isIORank) { return; } - - MessageBufferType buffer; - this->pack(0, buffer); - - // pass a dummy_link to satisfy virtual class - const int dummyLink = -1; - this->unpack(dummyLink, buffer); - } - - void pack(int link, MessageBufferType& buffer) - { - if (link != 0) - throw std::logic_error("link in method pack is not 0 as expected"); - for (int i = 0; i < 3; ++i) { - const auto& name = localFlows_[i].first; - buffer.write(name); - unsigned int size = localFlows_[i].second.first.size(); - buffer.write(size); - for (unsigned int j = 0; j < size; ++j) { - const auto& nncIdx = localFlows_[i].second.first[j]; - buffer.write(nncIdx); - const auto& flows = localFlows_[i].second.second[j]; - buffer.write(flows); - } - } - } - - void unpack(int /*link*/, MessageBufferType& buffer) - { - for (int i = 0; i < 3; ++i) { - std::string name; - buffer.read(name); - globalFlows_[i].first = name; - unsigned int size = 0; - buffer.read(size); - for (unsigned int j = 0; j < size; ++j) { - int nncIdx; - double data; - buffer.read(nncIdx); - buffer.read(data); - if (nncIdx < 0) - continue; - // This array is initialized in the collect(...) method below - globalFlows_[i].second.second[nncIdx] = data; - } - } - } -}; - -template -CollectDataToIORank:: -CollectDataToIORank(const Grid& grid, const EquilGrid* equilGrid, - const GridView& localGridView, - const Dune::CartesianIndexMapper& cartMapper, - const Dune::CartesianIndexMapper* equilCartMapper, - const std::set& fipRegionsInterregFlow) - : toIORankComm_(grid.comm()) - , globalInterRegFlows_(EclInterRegFlowMap::createMapFromNames(toVector(fipRegionsInterregFlow))) -{ - // index maps only have to be build when reordering is needed - if (!needsReordering && !isParallel()) - return; - - const CollectiveCommunication& comm = grid.comm(); - - { - std::set send, recv; - using EquilGridView = typename EquilGrid::LeafGridView; - typename std::is_same::type isSameGrid; - - typedef Dune::MultipleCodimMultipleGeomTypeMapper ElementMapper; - ElementMapper elemMapper(localGridView, Dune::mcmgElementLayout()); - sortedCartesianIdx_.reserve(localGridView.size(0)); - - for (const auto& elem : elements(localGridView)) - { - auto idx = elemMapper.index(elem); - sortedCartesianIdx_.push_back(cartMapper.cartesianIndex(idx)); - } - - std::sort(sortedCartesianIdx_.begin(), sortedCartesianIdx_.end()); - localIdxToGlobalIdx_.resize(localGridView.size(0), -1); - - // the I/O rank receives from all other ranks - if (isIORank()) { - // We need a mapping from local to global grid, here we - // use equilGrid which represents a view on the global grid - // reserve memory - const size_t globalSize = equilGrid->leafGridView().size(0); - globalCartesianIndex_.resize(globalSize, -1); - const EquilGridView equilGridView = equilGrid->leafGridView(); - - using EquilElementMapper = Dune::MultipleCodimMultipleGeomTypeMapper; - EquilElementMapper equilElemMapper(equilGridView, Dune::mcmgElementLayout()); - - // Scatter the global index to local index for lookup during restart - if constexpr (isSameGrid) { - ElementIndexScatterHandle handle(equilElemMapper, elemMapper, localIdxToGlobalIdx_); - grid.scatterData(handle); - } - - // loop over all elements (global grid) and store Cartesian index - for (const auto& elem : elements(equilGrid->leafGridView())) { - int elemIdx = equilElemMapper.index(elem); - int cartElemIdx = equilCartMapper->cartesianIndex(elemIdx); - globalCartesianIndex_[elemIdx] = cartElemIdx; - } - - for (int i = 0; i < comm.size(); ++i) { - if (i != ioRank) - recv.insert(i); - } - } - else - { - // all other simply send to the I/O rank - send.insert(ioRank); - - // Scatter the global index to local index for lookup during restart - // This is a bit hacky since the type differs from the iorank. - // But should work since we only receive, i.e. use the second parameter. - if constexpr (isSameGrid) { - ElementIndexScatterHandle handle(elemMapper, elemMapper, localIdxToGlobalIdx_); - grid.scatterData(handle); - } - } - - // Sync the global element indices - if constexpr (isSameGrid) { - ElementIndexHandle handle(elemMapper, localIdxToGlobalIdx_); - grid.communicate(handle, Dune::InteriorBorder_All_Interface, - Dune::ForwardCommunication); - } - - localIndexMap_.clear(); - const size_t gridSize = grid.size(0); - localIndexMap_.reserve(gridSize); - - // store the local Cartesian index - IndexMapType distributedCartesianIndex; - distributedCartesianIndex.resize(gridSize, -1); - - // A mapping for the whole grid (including the ghosts) is needed for restarts - for (const auto& elem : elements(localGridView)) { - int elemIdx = elemMapper.index(elem); - distributedCartesianIndex[elemIdx] = cartMapper.cartesianIndex(elemIdx); - - // only store interior element for collection - //assert(element.partitionType() == Dune::InteriorEntity); - - localIndexMap_.push_back(elemIdx); - } - - // insert send and recv linkage to communicator - toIORankComm_.insertRequest(send, recv); - - // need an index map for each rank - indexMaps_.clear(); - indexMaps_.resize(comm.size()); - - // distribute global id's to io rank for later association of dof's - DistributeIndexMapping distIndexMapping(globalCartesianIndex_, - distributedCartesianIndex, - localIndexMap_, - indexMaps_, - globalRanks_, - recv, - isIORank()); - toIORankComm_.exchange(distIndexMapping); - } -} - -template -void CollectDataToIORank:: -collect(const data::Solution& localCellData, - const std::map, double>& localBlockData, - const data::Wells& localWellData, - const data::WellBlockAveragePressures& localWBPData, - const data::GroupAndNetworkValues& localGroupAndNetworkData, - const data::Aquifers& localAquiferData, - const WellTestState& localWellTestState, - const EclInterRegFlowMap& localInterRegFlows, - const std::array, std::vector>>, 3>& localFlowsn, - const std::array, std::vector>>, 3>& localFloresn) -{ - globalCellData_ = {}; - globalBlockData_.clear(); - globalWellData_.clear(); - globalWBPData_.values.clear(); - globalGroupAndNetworkData_.clear(); - globalAquiferData_.clear(); - globalWellTestState_.clear(); - this->globalInterRegFlows_.clear(); - globalFlowsn_ = {}; - globalFloresn_ = {}; - - // index maps only have to be build when reordering is needed - if(!needsReordering && !isParallel()) - return; - - // this also linearises the local buffers on ioRank - PackUnPackCellData packUnpackCellData { - localCellData, - this->globalCellData_, - this->localIndexMap_, - this->indexMaps_, - this->numCells(), - this->isIORank() - }; - - if (! isParallel()) { - // no need to collect anything. - return; - } - - // Set the right sizes for Flowsn and Floresn - for (int i = 0; i < 3; ++i) { - unsigned int sizeFlr = localFloresn[i].second.first.size(); - globalFloresn_[i].second.first.resize(sizeFlr, 0); - globalFloresn_[i].second.second.resize(sizeFlr, 0.0); - unsigned int sizeFlo = localFlowsn[i].second.first.size(); - globalFlowsn_[i].second.first.resize(sizeFlo, 0); - globalFlowsn_[i].second.second.resize(sizeFlo, 0.0); - } - - PackUnPackWellData packUnpackWellData { - localWellData, - this->globalWellData_, - this->isIORank() - }; - - PackUnPackGroupAndNetworkValues packUnpackGroupAndNetworkData { - localGroupAndNetworkData, - this->globalGroupAndNetworkData_, - this->isIORank() - }; - - PackUnPackBlockData packUnpackBlockData { - localBlockData, - this->globalBlockData_, - this->isIORank() - }; - - PackUnPackWBPData packUnpackWBPData { - localWBPData, - this->globalWBPData_, - this->isIORank() - }; - - PackUnPackAquiferData packUnpackAquiferData { - localAquiferData, - this->globalAquiferData_, - this->isIORank() - }; - - PackUnPackWellTestState packUnpackWellTestState { - localWellTestState, - this->globalWellTestState_, - this->isIORank() - }; - - PackUnpackInterRegFlows packUnpackInterRegFlows { - localInterRegFlows, - this->globalInterRegFlows_, - this->isIORank() - }; - - PackUnpackFlows packUnpackFlowsn { - localFlowsn, - this->globalFlowsn_, - this->isIORank() - }; - - PackUnpackFlows packUnpackFloresn { - localFloresn, - this->globalFloresn_, - this->isIORank() - }; - - toIORankComm_.exchange(packUnpackCellData); - toIORankComm_.exchange(packUnpackWellData); - toIORankComm_.exchange(packUnpackGroupAndNetworkData); - toIORankComm_.exchange(packUnpackBlockData); - toIORankComm_.exchange(packUnpackWBPData); - toIORankComm_.exchange(packUnpackAquiferData); - toIORankComm_.exchange(packUnpackWellTestState); - toIORankComm_.exchange(packUnpackInterRegFlows); - toIORankComm_.exchange(packUnpackFlowsn); - toIORankComm_.exchange(packUnpackFloresn); - -#ifndef NDEBUG - // make sure every process is on the same page - toIORankComm_.barrier(); -#endif -} - -template -int CollectDataToIORank:: -localIdxToGlobalIdx(unsigned localIdx) const -{ - if (!isParallel()) { - return localIdx; - } - - if (this->localIdxToGlobalIdx_.empty()) { - throw std::logic_error("index map is not created on this rank"); - } - - if (localIdx > this->localIdxToGlobalIdx_.size()) { - throw std::logic_error("local index is outside map range"); - } - - return this->localIdxToGlobalIdx_[localIdx]; -} - -template -bool CollectDataToIORank:: -isCartIdxOnThisRank(int cartIdx) const -{ - if (! this->isParallel()) { - return true; - } - - assert (! needsReordering); - - auto candidate = std::lower_bound(this->sortedCartesianIdx_.begin(), - this->sortedCartesianIdx_.end(), - cartIdx); - - return (candidate != sortedCartesianIdx_.end()) - && (*candidate == cartIdx); -} - - - #if HAVE_DUNE_FEM template class CollectDataToIORank > >; - -#ifdef HAVE_DUNE_ALUGRID -#if HAVE_MPI - using ALUGrid3CN = Dune::ALUGrid<3, 3, Dune::cube, Dune::nonconforming, Dune::ALUGridMPIComm>; -#else - using ALUGrid3CN = Dune::ALUGrid<3, 3, Dune::cube, Dune::nonconforming, Dune::ALUGridNoComm>; -#endif //HAVE_MPI - -template class CollectDataToIORank>>>; - -template class CollectDataToIORank>>; - -#endif //HAVE_DUNE_ALUGRID - #else // ! HAVE_DUNE_FEM template class CollectDataToIORank>>; - -#ifdef HAVE_DUNE_ALUGRID -#if HAVE_MPI - using ALUGrid3CN = Dune::ALUGrid<3, 3, Dune::cube, Dune::nonconforming, Dune::ALUGridMPIComm>; -#else - using ALUGrid3CN = Dune::ALUGrid<3, 3, Dune::cube, Dune::nonconforming, Dune::ALUGridNoComm>; -#endif //HAVE_MPI -template class CollectDataToIORank>>; - -#endif //HAVE_DUNE_ALUGRID -#endif // ! HAVE_DUNE_FEM -template class CollectDataToIORank, - Dune::PolyhedralGrid<3,3,double>, - Dune::GridView>>; +#endif } // end namespace Opm diff --git a/ebos/collecttoiorank_impl.hh b/ebos/collecttoiorank_impl.hh new file mode 100644 index 000000000..f179b9c39 --- /dev/null +++ b/ebos/collecttoiorank_impl.hh @@ -0,0 +1,1128 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/* + 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 2 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 . + + Consult the COPYING file in the top-level source directory of this + module for the precise wording of the license and the list of + copyright holders. +*/ + +#ifndef EWOMS_COLLECT_TO_IO_RANK_IMPL_HH +#define EWOMS_COLLECT_TO_IO_RANK_IMPL_HH + +#include + +#include + +#include +#include +#include + +#include +#include +#include +#include +#include + +namespace { + std::vector toVector(const std::set& string_set) + { + return { string_set.begin(), string_set.end() }; + } +} + +namespace Opm { + +// global id +class GlobalCellIndex +{ + int globalId_; + int localIndex_; + bool isInterior_; + +public: + GlobalCellIndex() + : globalId_(-1) + , localIndex_(-1) + , isInterior_(true) + {} + void setGhost() + { isInterior_ = false; } + + void setId(int globalId) + { globalId_ = globalId; } + void setIndex(int localIndex) + { localIndex_ = localIndex; } + + int localIndex () const + { return localIndex_; } + int id () const + { return globalId_; } + bool isInterior() const + { return isInterior_; } +}; + +using IndexMapType = std::vector; +using IndexMapStorageType = std::vector; +using P2PCommunicatorType = Dune::Point2PointCommunicator; +using MessageBufferType = typename P2PCommunicatorType::MessageBufferType; +class DistributeIndexMapping : public P2PCommunicatorType::DataHandleInterface +{ +protected: + const std::vector& distributedGlobalIndex_; + IndexMapType& localIndexMap_; + IndexMapStorageType& indexMaps_; + std::map globalPosition_; + std::set& recv_; + std::vector& ranks_; + +public: + DistributeIndexMapping(const std::vector& globalIndex, + const std::vector& distributedGlobalIndex, + IndexMapType& localIndexMap, + IndexMapStorageType& indexMaps, + std::vector& ranks, + std::set& recv, + bool isIORank) + : distributedGlobalIndex_(distributedGlobalIndex) + , localIndexMap_(localIndexMap) + , indexMaps_(indexMaps) + , globalPosition_() + , recv_(recv) + , ranks_(ranks) + { + size_t size = globalIndex.size(); + // create mapping globalIndex --> localIndex + if ( isIORank ) // ioRank + for (size_t index = 0; index < size; ++index) + globalPosition_.insert(std::make_pair(globalIndex[index], index)); + + // we need to create a mapping from local to global + if (!indexMaps_.empty()) { + if (isIORank) + ranks_.resize(size, -1); + // for the ioRank create a localIndex to index in global state map + IndexMapType& indexMap = indexMaps_.back(); + size_t localSize = localIndexMap_.size(); + indexMap.resize(localSize); + for (size_t i=0; isecond; + // Using max should be backwards compatible + ranks_[entry] = std::max(ranks_[entry], rank); + } + if (rankIt != recv_.end()) + ++rankIt; + } + #ifndef NDEBUG + for (const auto& rank: ranks_) + assert(rank>=0); + #endif + } + } + + void pack(int link, MessageBufferType& buffer) + { + // we should only get one link + if (link != 0) + throw std::logic_error("link in method pack is not 0 as execpted"); + + // pack all interior global cell id's + int size = localIndexMap_.size(); + buffer.write(size); + + for (int index = 0; index < size; ++index) { + int globalIdx = distributedGlobalIndex_[localIndexMap_[index]]; + buffer.write(globalIdx); + } + } + + void unpack(int link, MessageBufferType& buffer) + { + // get index map for current link + IndexMapType& indexMap = indexMaps_[link]; + assert(!globalPosition_.empty()); + + // unpack all interior global cell id's + int numCells = 0; + buffer.read(numCells); + indexMap.resize(numCells); + for (int index = 0; index < numCells; ++index) { + buffer.read(indexMap[index]); + } + } +}; + + +/// \brief Communication handle to scatter the global index +template +class ElementIndexScatterHandle +{ +public: + ElementIndexScatterHandle(const EquilMapper& sendMapper, const Mapper& recvMapper, std::vector& elementIndices) + : sendMapper_(sendMapper), recvMapper_(recvMapper), elementIndices_(elementIndices) + {} + using DataType = int; + bool fixedSize(int /*dim*/, int /*codim*/) + { + return true; + } + + template + std::size_t size(const T&) + { + return 1; + } + template + void gather(B& buffer, const T& t) + { + buffer.write(sendMapper_.index(t)); + } + template + void scatter(B& buffer, const T& t, std::size_t) + { + buffer.read(elementIndices_[recvMapper_.index(t)]); + } + + bool contains(int dim, int codim) + { + return dim==3 && codim==0; + } +private: + const EquilMapper& sendMapper_; + const Mapper& recvMapper_; + std::vector& elementIndices_; +}; + +/// \brief Communication handle to scatter the global index +template +class ElementIndexHandle +{ +public: + ElementIndexHandle(const Mapper& mapper, std::vector& elementIndices) + : mapper_(mapper), elementIndices_(elementIndices) + {} + using DataType = int; + bool fixedSize(int /*dim*/, int /*codim*/) + { + return true; + } + + template + std::size_t size(const T&) + { + return 1; + } + template + void gather(B& buffer, const T& t) + { + buffer.write(elementIndices_[mapper_.index(t)]); + } + template + void scatter(B& buffer, const T& t, std::size_t) + { + buffer.read(elementIndices_[mapper_.index(t)]); + } + + bool contains(int dim, int codim) + { + return dim==3 && codim==0; + } +private: + const Mapper& mapper_; + std::vector& elementIndices_; +}; + +class PackUnPackCellData : public P2PCommunicatorType::DataHandleInterface +{ + const data::Solution& localCellData_; + data::Solution& globalCellData_; + + const IndexMapType& localIndexMap_; + const IndexMapStorageType& indexMaps_; + +public: + PackUnPackCellData(const data::Solution& localCellData, + data::Solution& globalCellData, + const IndexMapType& localIndexMap, + const IndexMapStorageType& indexMaps, + size_t globalSize, + bool isIORank) + : localCellData_(localCellData) + , globalCellData_(globalCellData) + , localIndexMap_(localIndexMap) + , indexMaps_(indexMaps) + { + if (isIORank) { + // add missing data to global cell data + for (const auto& pair : localCellData_) { + const std::string& key = pair.first; + std::size_t containerSize = globalSize; + [[maybe_unused]] auto ret = globalCellData_.insert(key, pair.second.dim, + std::vector(containerSize), + pair.second.target); + assert(ret.second); + } + + MessageBufferType buffer; + pack(0, buffer); + + // the last index map is the local one + doUnpack(indexMaps.back(), buffer); + } + } + + // pack all data associated with link + void pack(int link, MessageBufferType& buffer) + { + // we should only get one link + if (link != 0) + throw std::logic_error("link in method pack is not 0 as expected"); + + // write all cell data registered in local state + for (const auto& pair : localCellData_) { + const auto& data = pair.second.data; + + // write all data from local data to buffer + write(buffer, localIndexMap_, data); + } + } + + void doUnpack(const IndexMapType& indexMap, MessageBufferType& buffer) + { + // we loop over the data as + // its order governs the order the data got received. + for (auto& pair : localCellData_) { + const std::string& key = pair.first; + auto& data = globalCellData_.data(key); + + //write all data from local cell data to buffer + read(buffer, indexMap, data); + } + } + + // unpack all data associated with link + void unpack(int link, MessageBufferType& buffer) + { doUnpack(indexMaps_[link], buffer); } + +protected: + template + void write(MessageBufferType& buffer, + const IndexMapType& localIndexMap, + const Vector& vector, + unsigned int offset = 0, + unsigned int stride = 1) const + { + unsigned int size = localIndexMap.size(); + buffer.write(size); + assert(vector.size() >= stride * size); + for (unsigned int i=0; i + void read(MessageBufferType& buffer, + const IndexMapType& indexMap, + Vector& vector, + unsigned int offset = 0, + unsigned int stride = 1) const + { + unsigned int size = 0; + buffer.read(size); + assert(size == indexMap.size()); + for (unsigned int i=0; ipack(0, buffer); + + // pass a dummy_link to satisfy virtual class + const int dummyLink = -1; + this->unpack(dummyLink, buffer); + } + + // pack all data associated with link + void pack(int link, MessageBufferType& buffer) + { + // we should only get one link + if (link != 0) { + throw std::logic_error { + "link in method pack is not 0 as expected" + }; + } + + // write all group and network (node/branch) data + this->localGroupAndNetworkData_.write(buffer); + } + + // unpack all data associated with link + void unpack(int /*link*/, MessageBufferType& buffer) + { this->globalGroupAndNetworkData_.read(buffer); } +}; + +class PackUnPackBlockData : public P2PCommunicatorType::DataHandleInterface +{ + const std::map, double>& localBlockData_; + std::map, double>& globalBlockValues_; + +public: + PackUnPackBlockData(const std::map, double>& localBlockData, + std::map, double>& globalBlockValues, + bool isIORank) + : localBlockData_(localBlockData) + , globalBlockValues_(globalBlockValues) + { + if (isIORank) { + MessageBufferType buffer; + pack(0, buffer); + + // pass a dummyLink to satisfy virtual class + int dummyLink = -1; + unpack(dummyLink, buffer); + } + } + + // pack all data associated with link + void pack(int link, MessageBufferType& buffer) + { + // we should only get one link + if (link != 0) + throw std::logic_error("link in method pack is not 0 as expected"); + + // write all block data + unsigned int size = localBlockData_.size(); + buffer.write(size); + for (const auto& map : localBlockData_) { + buffer.write(map.first.first); + buffer.write(map.first.second); + buffer.write(map.second); + } + } + + // unpack all data associated with link + void unpack(int /*link*/, MessageBufferType& buffer) + { + // read all block data + unsigned int size = 0; + buffer.read(size); + for (size_t i = 0; i < size; ++i) { + std::string name; + int idx; + double data; + buffer.read(name); + buffer.read(idx); + buffer.read(data); + globalBlockValues_[std::make_pair(name, idx)] = data; + } + } +}; + +class PackUnPackWBPData : public P2PCommunicatorType::DataHandleInterface +{ + const data::WellBlockAveragePressures& localWBPData_; + data::WellBlockAveragePressures& globalWBPValues_; + +public: + PackUnPackWBPData(const data::WellBlockAveragePressures& localWBPData, + data::WellBlockAveragePressures& globalWBPValues, + const bool isIORank) + : localWBPData_ (localWBPData) + , globalWBPValues_(globalWBPValues) + { + if (! isIORank) { + return; + } + + MessageBufferType buffer; + pack(0, buffer); + + // Pass a dummy link to satisfy base class API requirement + const int dummyLink = -1; + unpack(dummyLink, buffer); + } + + // Pack all data associated with link + void pack(int link, MessageBufferType& buffer) + { + // We should only get one link + if (link != 0) { + throw std::logic_error { + "link (" + std::to_string(link) + + ") in method pack() is not 0 as expected" + }; + } + + // Write all local, per-well, WBP data + this->localWBPData_.write(buffer); + } + + // Unpack all data associated with link + void unpack([[maybe_unused]] const int link, + MessageBufferType& buffer) + { + this->globalWBPValues_.read(buffer); + } +}; + +class PackUnPackWellTestState : public P2PCommunicatorType::DataHandleInterface +{ +public: + PackUnPackWellTestState(const WellTestState& localWTestState, + WellTestState& globalWTestState, + bool isIORank) + : local_(localWTestState) + , global_(globalWTestState) + { + if (isIORank) { + MessageBufferType buffer; + pack(0, buffer); + + // pass a dummyLink to satisfy virtual class + int dummyLink = -1; + unpack(dummyLink, buffer); + } + } + + void pack(int link, MessageBufferType& buffer) { + if (link != 0) + throw std::logic_error("link in method pack is not 0 as expected"); + this->local_.pack(buffer); + } + + void unpack(int, MessageBufferType& buffer) { + this->global_.unpack(buffer); + } + +private: + const WellTestState& local_; + WellTestState& global_; +}; + +class PackUnPackAquiferData : public P2PCommunicatorType::DataHandleInterface +{ + const data::Aquifers& localAquiferData_; + data::Aquifers& globalAquiferData_; + +public: + PackUnPackAquiferData(const data::Aquifers& localAquiferData, + data::Aquifers& globalAquiferData, + bool isIORank) + : localAquiferData_(localAquiferData) + , globalAquiferData_(globalAquiferData) + { + if (isIORank) { + MessageBufferType buffer; + pack(0, buffer); + + // pass a dummyLink to satisfy virtual class + int dummyLink = -1; + unpack(dummyLink, buffer); + } + } + + // pack all data associated with link + void pack(int link, MessageBufferType& buffer) + { + // we should only get one link + if (link != 0) + throw std::logic_error("link in method pack is not 0 as expected"); + + int size = localAquiferData_.size(); + buffer.write(size); + for (const auto& [key, data] : localAquiferData_) { + buffer.write(key); + data.write(buffer); + } + } + + // unpack all data associated with link + void unpack(int /*link*/, MessageBufferType& buffer) + { + int size; + buffer.read(size); + for (int i = 0; i < size; ++i) { + int key; + buffer.read(key); + data::AquiferData data; + data.read(buffer); + + auto& aq = this->globalAquiferData_[key]; + + this->unpackCommon(data, aq); + + if (auto const* aquFet = data.typeData.get(); + aquFet != nullptr) + { + this->unpackAquFet(*aquFet, aq); + } + else if (auto const* aquCT = data.typeData.get(); + aquCT != nullptr) + { + this->unpackCarterTracy(*aquCT, aq); + } + else if (auto const* aquNum = data.typeData.get(); + aquNum != nullptr) + { + this->unpackNumericAquifer(*aquNum, aq); + } + } + } + +private: + void unpackCommon(const data::AquiferData& data, data::AquiferData& aq) + { + aq.aquiferID = std::max(aq.aquiferID, data.aquiferID); + aq.pressure = std::max(aq.pressure, data.pressure); + aq.initPressure = std::max(aq.initPressure, data.initPressure); + aq.datumDepth = std::max(aq.datumDepth, data.datumDepth); + aq.fluxRate += data.fluxRate; + aq.volume += data.volume; + } + + void unpackAquFet(const data::FetkovichData& aquFet, data::AquiferData& aq) + { + if (! aq.typeData.is()) { + auto* tData = aq.typeData.create(); + *tData = aquFet; + } + else { + const auto& src = aquFet; + auto& dst = *aq.typeData.getMutable(); + + dst.initVolume = std::max(dst.initVolume , src.initVolume); + dst.prodIndex = std::max(dst.prodIndex , src.prodIndex); + dst.timeConstant = std::max(dst.timeConstant, src.timeConstant); + } + } + + void unpackCarterTracy(const data::CarterTracyData& aquCT, data::AquiferData& aq) + { + if (! aq.typeData.is()) { + auto* tData = aq.typeData.create(); + *tData = aquCT; + } + else { + const auto& src = aquCT; + auto& dst = *aq.typeData.getMutable(); + + dst.timeConstant = std::max(dst.timeConstant , src.timeConstant); + dst.influxConstant = std::max(dst.influxConstant, src.influxConstant); + dst.waterDensity = std::max(dst.waterDensity , src.waterDensity); + dst.waterViscosity = std::max(dst.waterViscosity, src.waterViscosity); + + dst.dimensionless_time = std::max(dst.dimensionless_time , src.dimensionless_time); + dst.dimensionless_pressure = std::max(dst.dimensionless_pressure, src.dimensionless_pressure); + } + } + + void unpackNumericAquifer(const data::NumericAquiferData& aquNum, data::AquiferData& aq) + { + if (! aq.typeData.is()) { + auto* tData = aq.typeData.create(); + *tData = aquNum; + } + else { + const auto& src = aquNum; + auto& dst = *aq.typeData.getMutable(); + + std::transform(src.initPressure.begin(), + src.initPressure.end(), + dst.initPressure.begin(), + dst.initPressure.begin(), + [](const double p0_1, const double p0_2) + { + return std::max(p0_1, p0_2); + }); + } + } +}; + +class PackUnpackInterRegFlows : public P2PCommunicatorType::DataHandleInterface +{ + const EclInterRegFlowMap& localInterRegFlows_; + EclInterRegFlowMap& globalInterRegFlows_; + +public: + PackUnpackInterRegFlows(const EclInterRegFlowMap& localInterRegFlows, + EclInterRegFlowMap& globalInterRegFlows, + const bool isIORank) + : localInterRegFlows_(localInterRegFlows) + , globalInterRegFlows_(globalInterRegFlows) + { + if (! isIORank) { return; } + + MessageBufferType buffer; + this->pack(0, buffer); + + // pass a dummy_link to satisfy virtual class + const int dummyLink = -1; + this->unpack(dummyLink, buffer); + } + + // pack all data associated with link + void pack(int link, MessageBufferType& buffer) + { + // we should only get one link + if (link != 0) { + throw std::logic_error { + "link in method pack is not 0 as expected" + }; + } + + // write all inter-region flow data + this->localInterRegFlows_.write(buffer); + } + + // unpack all data associated with link + void unpack(int /*link*/, MessageBufferType& buffer) + { this->globalInterRegFlows_.read(buffer); } +}; + +class PackUnpackFlows : public P2PCommunicatorType::DataHandleInterface +{ + const std::array, std::vector>>, 3>& localFlows_; + std::array, std::vector>>, 3>& globalFlows_; + +public: + PackUnpackFlows(const std::array, std::vector>>, 3> & localFlows, + std::array, std::vector>>, 3>& globalFlows, + const bool isIORank) + : localFlows_(localFlows) + , globalFlows_(globalFlows) + { + if (! isIORank) { return; } + + MessageBufferType buffer; + this->pack(0, buffer); + + // pass a dummy_link to satisfy virtual class + const int dummyLink = -1; + this->unpack(dummyLink, buffer); + } + + void pack(int link, MessageBufferType& buffer) + { + if (link != 0) + throw std::logic_error("link in method pack is not 0 as expected"); + for (int i = 0; i < 3; ++i) { + const auto& name = localFlows_[i].first; + buffer.write(name); + unsigned int size = localFlows_[i].second.first.size(); + buffer.write(size); + for (unsigned int j = 0; j < size; ++j) { + const auto& nncIdx = localFlows_[i].second.first[j]; + buffer.write(nncIdx); + const auto& flows = localFlows_[i].second.second[j]; + buffer.write(flows); + } + } + } + + void unpack(int /*link*/, MessageBufferType& buffer) + { + for (int i = 0; i < 3; ++i) { + std::string name; + buffer.read(name); + globalFlows_[i].first = name; + unsigned int size = 0; + buffer.read(size); + for (unsigned int j = 0; j < size; ++j) { + int nncIdx; + double data; + buffer.read(nncIdx); + buffer.read(data); + if (nncIdx < 0) + continue; + // This array is initialized in the collect(...) method below + globalFlows_[i].second.second[nncIdx] = data; + } + } + } +}; + +template +CollectDataToIORank:: +CollectDataToIORank(const Grid& grid, const EquilGrid* equilGrid, + const GridView& localGridView, + const Dune::CartesianIndexMapper& cartMapper, + const Dune::CartesianIndexMapper* equilCartMapper, + const std::set& fipRegionsInterregFlow) + : toIORankComm_(grid.comm()) + , globalInterRegFlows_(EclInterRegFlowMap::createMapFromNames(toVector(fipRegionsInterregFlow))) +{ + // index maps only have to be build when reordering is needed + if (!needsReordering && !isParallel()) + return; + + const CollectiveCommunication& comm = grid.comm(); + + { + std::set send, recv; + using EquilGridView = typename EquilGrid::LeafGridView; + typename std::is_same::type isSameGrid; + + typedef Dune::MultipleCodimMultipleGeomTypeMapper ElementMapper; + ElementMapper elemMapper(localGridView, Dune::mcmgElementLayout()); + sortedCartesianIdx_.reserve(localGridView.size(0)); + + for (const auto& elem : elements(localGridView)) + { + auto idx = elemMapper.index(elem); + sortedCartesianIdx_.push_back(cartMapper.cartesianIndex(idx)); + } + + std::sort(sortedCartesianIdx_.begin(), sortedCartesianIdx_.end()); + localIdxToGlobalIdx_.resize(localGridView.size(0), -1); + + // the I/O rank receives from all other ranks + if (isIORank()) { + // We need a mapping from local to global grid, here we + // use equilGrid which represents a view on the global grid + // reserve memory + const size_t globalSize = equilGrid->leafGridView().size(0); + globalCartesianIndex_.resize(globalSize, -1); + const EquilGridView equilGridView = equilGrid->leafGridView(); + + using EquilElementMapper = Dune::MultipleCodimMultipleGeomTypeMapper; + EquilElementMapper equilElemMapper(equilGridView, Dune::mcmgElementLayout()); + + // Scatter the global index to local index for lookup during restart + if constexpr (isSameGrid) { + ElementIndexScatterHandle handle(equilElemMapper, elemMapper, localIdxToGlobalIdx_); + grid.scatterData(handle); + } + + // loop over all elements (global grid) and store Cartesian index + for (const auto& elem : elements(equilGrid->leafGridView())) { + int elemIdx = equilElemMapper.index(elem); + int cartElemIdx = equilCartMapper->cartesianIndex(elemIdx); + globalCartesianIndex_[elemIdx] = cartElemIdx; + } + + for (int i = 0; i < comm.size(); ++i) { + if (i != ioRank) + recv.insert(i); + } + } + else + { + // all other simply send to the I/O rank + send.insert(ioRank); + + // Scatter the global index to local index for lookup during restart + // This is a bit hacky since the type differs from the iorank. + // But should work since we only receive, i.e. use the second parameter. + if constexpr (isSameGrid) { + ElementIndexScatterHandle handle(elemMapper, elemMapper, localIdxToGlobalIdx_); + grid.scatterData(handle); + } + } + + // Sync the global element indices + if constexpr (isSameGrid) { + ElementIndexHandle handle(elemMapper, localIdxToGlobalIdx_); + grid.communicate(handle, Dune::InteriorBorder_All_Interface, + Dune::ForwardCommunication); + } + + localIndexMap_.clear(); + const size_t gridSize = grid.size(0); + localIndexMap_.reserve(gridSize); + + // store the local Cartesian index + IndexMapType distributedCartesianIndex; + distributedCartesianIndex.resize(gridSize, -1); + + // A mapping for the whole grid (including the ghosts) is needed for restarts + for (const auto& elem : elements(localGridView)) { + int elemIdx = elemMapper.index(elem); + distributedCartesianIndex[elemIdx] = cartMapper.cartesianIndex(elemIdx); + + // only store interior element for collection + //assert(element.partitionType() == Dune::InteriorEntity); + + localIndexMap_.push_back(elemIdx); + } + + // insert send and recv linkage to communicator + toIORankComm_.insertRequest(send, recv); + + // need an index map for each rank + indexMaps_.clear(); + indexMaps_.resize(comm.size()); + + // distribute global id's to io rank for later association of dof's + DistributeIndexMapping distIndexMapping(globalCartesianIndex_, + distributedCartesianIndex, + localIndexMap_, + indexMaps_, + globalRanks_, + recv, + isIORank()); + toIORankComm_.exchange(distIndexMapping); + } +} + +template +void CollectDataToIORank:: +collect(const data::Solution& localCellData, + const std::map, double>& localBlockData, + const data::Wells& localWellData, + const data::WellBlockAveragePressures& localWBPData, + const data::GroupAndNetworkValues& localGroupAndNetworkData, + const data::Aquifers& localAquiferData, + const WellTestState& localWellTestState, + const EclInterRegFlowMap& localInterRegFlows, + const std::array, std::vector>>, 3>& localFlowsn, + const std::array, std::vector>>, 3>& localFloresn) +{ + globalCellData_ = {}; + globalBlockData_.clear(); + globalWellData_.clear(); + globalWBPData_.values.clear(); + globalGroupAndNetworkData_.clear(); + globalAquiferData_.clear(); + globalWellTestState_.clear(); + this->globalInterRegFlows_.clear(); + globalFlowsn_ = {}; + globalFloresn_ = {}; + + // index maps only have to be build when reordering is needed + if(!needsReordering && !isParallel()) + return; + + // this also linearises the local buffers on ioRank + PackUnPackCellData packUnpackCellData { + localCellData, + this->globalCellData_, + this->localIndexMap_, + this->indexMaps_, + this->numCells(), + this->isIORank() + }; + + if (! isParallel()) { + // no need to collect anything. + return; + } + + // Set the right sizes for Flowsn and Floresn + for (int i = 0; i < 3; ++i) { + unsigned int sizeFlr = localFloresn[i].second.first.size(); + globalFloresn_[i].second.first.resize(sizeFlr, 0); + globalFloresn_[i].second.second.resize(sizeFlr, 0.0); + unsigned int sizeFlo = localFlowsn[i].second.first.size(); + globalFlowsn_[i].second.first.resize(sizeFlo, 0); + globalFlowsn_[i].second.second.resize(sizeFlo, 0.0); + } + + PackUnPackWellData packUnpackWellData { + localWellData, + this->globalWellData_, + this->isIORank() + }; + + PackUnPackGroupAndNetworkValues packUnpackGroupAndNetworkData { + localGroupAndNetworkData, + this->globalGroupAndNetworkData_, + this->isIORank() + }; + + PackUnPackBlockData packUnpackBlockData { + localBlockData, + this->globalBlockData_, + this->isIORank() + }; + + PackUnPackWBPData packUnpackWBPData { + localWBPData, + this->globalWBPData_, + this->isIORank() + }; + + PackUnPackAquiferData packUnpackAquiferData { + localAquiferData, + this->globalAquiferData_, + this->isIORank() + }; + + PackUnPackWellTestState packUnpackWellTestState { + localWellTestState, + this->globalWellTestState_, + this->isIORank() + }; + + PackUnpackInterRegFlows packUnpackInterRegFlows { + localInterRegFlows, + this->globalInterRegFlows_, + this->isIORank() + }; + + PackUnpackFlows packUnpackFlowsn { + localFlowsn, + this->globalFlowsn_, + this->isIORank() + }; + + PackUnpackFlows packUnpackFloresn { + localFloresn, + this->globalFloresn_, + this->isIORank() + }; + + toIORankComm_.exchange(packUnpackCellData); + toIORankComm_.exchange(packUnpackWellData); + toIORankComm_.exchange(packUnpackGroupAndNetworkData); + toIORankComm_.exchange(packUnpackBlockData); + toIORankComm_.exchange(packUnpackWBPData); + toIORankComm_.exchange(packUnpackAquiferData); + toIORankComm_.exchange(packUnpackWellTestState); + toIORankComm_.exchange(packUnpackInterRegFlows); + toIORankComm_.exchange(packUnpackFlowsn); + toIORankComm_.exchange(packUnpackFloresn); + +#ifndef NDEBUG + // make sure every process is on the same page + toIORankComm_.barrier(); +#endif +} + +template +int CollectDataToIORank:: +localIdxToGlobalIdx(unsigned localIdx) const +{ + if (!isParallel()) { + return localIdx; + } + + if (this->localIdxToGlobalIdx_.empty()) { + throw std::logic_error("index map is not created on this rank"); + } + + if (localIdx > this->localIdxToGlobalIdx_.size()) { + throw std::logic_error("local index is outside map range"); + } + + return this->localIdxToGlobalIdx_[localIdx]; +} + +template +bool CollectDataToIORank:: +isCartIdxOnThisRank(int cartIdx) const +{ + if (! this->isParallel()) { + return true; + } + + assert (! needsReordering); + + auto candidate = std::lower_bound(this->sortedCartesianIdx_.begin(), + this->sortedCartesianIdx_.end(), + cartIdx); + + return (candidate != sortedCartesianIdx_.end()) + && (*candidate == cartIdx); +} + +} // end namespace Opm +#endif diff --git a/ebos/eclgenericproblem.cc b/ebos/eclgenericproblem.cc index a4944c9a4..ce133dadb 100644 --- a/ebos/eclgenericproblem.cc +++ b/ebos/eclgenericproblem.cc @@ -20,30 +20,17 @@ module for the precise wording of the license and the list of copyright holders. */ - #include + #include +#include -#include - -#include -#include -#include - -#include - -#include +#include +#include #include -#include -#include - -#if HAVE_DUNE_ALUGRID -#include -#include -#include "alucartesianindexmapper.hh" -#endif // HAVE_DUNE_ALUGRID +#include #if HAVE_DUNE_FEM #include @@ -51,843 +38,23 @@ #include #endif // HAVE_DUNE_FEM -#include - -#include -#include -#include -#include - namespace Opm { -int eclPositionalParameter(Dune::ParameterTree& tree, - std::set& seenParams, - std::string& errorMsg, - const char** argv, - int paramIdx) -{ - std::string param = argv[paramIdx]; - size_t i = param.find('='); - if (i != std::string::npos) { - std::string oldParamName = param.substr(0, i); - std::string oldParamValue = param.substr(i+1); - std::string newParamName = "--" + oldParamName; - std::replace(newParamName.begin(), - newParamName.end(), '_' , '-'); - errorMsg = - "The old syntax to specify parameters on the command line is no longer supported: " - "Try replacing '" + oldParamName + "=" + oldParamValue + "' with "+ - "'" + newParamName + "=" + oldParamValue + "'!"; - return 0; - } - - if (seenParams.count("EclDeckFileName") > 0) { - errorMsg = - "Parameter 'EclDeckFileName' specified multiple times" - " as a command line parameter"; - return 0; - } - - tree["EclDeckFileName"] = argv[paramIdx]; - seenParams.insert("EclDeckFileName"); - return 1; -} - -template -EclGenericProblem:: -EclGenericProblem(const EclipseState& eclState, - const Schedule& schedule, - const GridView& gridView) - : eclState_(eclState) - , schedule_(schedule) - , gridView_(gridView) -{ -} - -template -EclGenericProblem -EclGenericProblem:: -serializationTestObject(const EclipseState& eclState, - const Schedule& schedule, - const GridView& gridView) -{ - EclGenericProblem result(eclState, schedule, gridView); - result.maxOilSaturation_ = {1.0, 2.0}; - result.maxWaterSaturation_ = {6.0}; - result.minOilPressure_ = {7.0, 8.0, 9.0, 10.0}; - result.overburdenPressure_ = {11.0}; - result.solventSaturation_ = {15.0}; - result.polymer_ = PolymerSolutionContainer::serializationTestObject(); - result.micp_ = MICPSolutionContainer::serializationTestObject(); - result.lastRv_ = {21.0}; - result.maxDRv_ = {22.0, 23.0}; - result.convectiveDrs_ = {24.0, 25.0, 26.0}; - result.lastRs_ = {27.0}; - result.maxDRs_ = {28.0}; - result.dRsDtOnlyFreeGas_ = {false, true}; - - return result; -} - -template -std::string -EclGenericProblem:: -helpPreamble(int, - const char **argv) -{ - std::string desc = EclGenericProblem::briefDescription(); - if (!desc.empty()) - desc = desc + "\n"; - - return - "Usage: "+std::string(argv[0]) + " [OPTIONS] [ECL_DECK_FILENAME]\n" - + desc; -} - -template -std::string -EclGenericProblem:: -briefDescription() -{ - if (briefDescription_.empty()) - return - "The Ecl-deck Black-Oil reservoir Simulator (ebos); a hydrocarbon " - "reservoir simulation program that processes ECL-formatted input " - "files that is part of the Open Porous Media project " - "(https://opm-project.org).\n" - "\n" - "THE GOAL OF THE `ebos` SIMULATOR IS TO CATER FOR THE NEEDS OF " - "DEVELOPMENT AND RESEARCH. No guarantees are made for production use!"; - else - return briefDescription_; -} - -template -void EclGenericProblem:: -readRockParameters_(const std::vector& cellCenterDepths) -{ - const auto& rock_config = eclState_.getSimulationConfig().rock_config(); - - // read the rock compressibility parameters - { - const auto& comp = rock_config.comp(); - rockParams_.clear(); - for (const auto& c : comp) - rockParams_.push_back( { c.pref, c.compressibility } ); - } - - // read the parameters for water-induced rock compaction - readRockCompactionParameters_(); - - unsigned numElem = gridView_.size(0); - if (eclState_.fieldProps().has_int(rock_config.rocknum_property())) { - rockTableIdx_.resize(numElem); - const auto& num = eclState_.fieldProps().get_int(rock_config.rocknum_property()); - for (size_t elemIdx = 0; elemIdx < numElem; ++ elemIdx) { - rockTableIdx_[elemIdx] = num[elemIdx] - 1; - } - } - - // Store overburden pressure pr element - const auto& overburdTables = eclState_.getTableManager().getOverburdTables(); - if (!overburdTables.empty()) { - overburdenPressure_.resize(numElem,0.0); - size_t numRocktabTables = rock_config.num_rock_tables(); - - if (overburdTables.size() != numRocktabTables) - throw std::runtime_error(std::to_string(numRocktabTables) +" OVERBURD tables is expected, but " + std::to_string(overburdTables.size()) +" is provided"); - - std::vector> overburdenTables(numRocktabTables); - for (size_t regionIdx = 0; regionIdx < numRocktabTables; ++regionIdx) { - const OverburdTable& overburdTable = overburdTables.template getTable(regionIdx); - overburdenTables[regionIdx].setXYContainers(overburdTable.getDepthColumn(),overburdTable.getOverburdenPressureColumn()); - } - - for (size_t elemIdx = 0; elemIdx < numElem; ++ elemIdx) { - unsigned tableIdx = 0; - if (!rockTableIdx_.empty()) { - tableIdx = rockTableIdx_[elemIdx]; - } - overburdenPressure_[elemIdx] = overburdenTables[tableIdx].eval(cellCenterDepths[elemIdx], /*extrapolation=*/true); - } - } -} - -template -void EclGenericProblem:: -readRockCompactionParameters_() -{ - const auto& rock_config = eclState_.getSimulationConfig().rock_config(); - - if (!rock_config.active()) - return; // deck does not enable rock compaction - - unsigned numElem = gridView_.size(0); - switch (rock_config.hysteresis_mode()) { - case RockConfig::Hysteresis::REVERS: - break; - case RockConfig::Hysteresis::IRREVERS: - // interpolate the porv volume multiplier using the minimum pressure in the cell - // i.e. don't allow re-inflation. - minOilPressure_.resize(numElem, 1e99); - break; - default: - throw std::runtime_error("Not support ROCKOMP hysteresis option "); - } - - size_t numRocktabTables = rock_config.num_rock_tables(); - bool waterCompaction = rock_config.water_compaction(); - - if (!waterCompaction) { - const auto& rocktabTables = eclState_.getTableManager().getRocktabTables(); - if (rocktabTables.size() != numRocktabTables) - throw std::runtime_error("ROCKCOMP is activated." + std::to_string(numRocktabTables) - +" ROCKTAB tables is expected, but " + std::to_string(rocktabTables.size()) +" is provided"); - - rockCompPoroMult_.resize(numRocktabTables); - rockCompTransMult_.resize(numRocktabTables); - for (size_t regionIdx = 0; regionIdx < numRocktabTables; ++regionIdx) { - const auto& rocktabTable = rocktabTables.template getTable(regionIdx); - const auto& pressureColumn = rocktabTable.getPressureColumn(); - const auto& poroColumn = rocktabTable.getPoreVolumeMultiplierColumn(); - const auto& transColumn = rocktabTable.getTransmissibilityMultiplierColumn(); - rockCompPoroMult_[regionIdx].setXYContainers(pressureColumn, poroColumn); - rockCompTransMult_[regionIdx].setXYContainers(pressureColumn, transColumn); - } - } else { - const auto& rock2dTables = eclState_.getTableManager().getRock2dTables(); - const auto& rock2dtrTables = eclState_.getTableManager().getRock2dtrTables(); - const auto& rockwnodTables = eclState_.getTableManager().getRockwnodTables(); - maxWaterSaturation_.resize(numElem, 0.0); - - if (rock2dTables.size() != numRocktabTables) - throw std::runtime_error("Water compation option is selected in ROCKCOMP." + std::to_string(numRocktabTables) - +" ROCK2D tables is expected, but " + std::to_string(rock2dTables.size()) +" is provided"); - - if (rockwnodTables.size() != numRocktabTables) - throw std::runtime_error("Water compation option is selected in ROCKCOMP." + std::to_string(numRocktabTables) - +" ROCKWNOD tables is expected, but " + std::to_string(rockwnodTables.size()) +" is provided"); - //TODO check size match - rockCompPoroMultWc_.resize(numRocktabTables, TabulatedTwoDFunction(TabulatedTwoDFunction::InterpolationPolicy::Vertical)); - for (size_t regionIdx = 0; regionIdx < numRocktabTables; ++regionIdx) { - const RockwnodTable& rockwnodTable = rockwnodTables.template getTable(regionIdx); - const auto& rock2dTable = rock2dTables[regionIdx]; - - if (rockwnodTable.getSaturationColumn().size() != rock2dTable.sizeMultValues()) - throw std::runtime_error("Number of entries in ROCKWNOD and ROCK2D needs to match."); - - for (size_t xIdx = 0; xIdx < rock2dTable.size(); ++xIdx) { - rockCompPoroMultWc_[regionIdx].appendXPos(rock2dTable.getPressureValue(xIdx)); - for (size_t yIdx = 0; yIdx < rockwnodTable.getSaturationColumn().size(); ++yIdx) - rockCompPoroMultWc_[regionIdx].appendSamplePoint(xIdx, - rockwnodTable.getSaturationColumn()[yIdx], - rock2dTable.getPvmultValue(xIdx, yIdx)); - } - } - - if (!rock2dtrTables.empty()) { - rockCompTransMultWc_.resize(numRocktabTables, TabulatedTwoDFunction(TabulatedTwoDFunction::InterpolationPolicy::Vertical)); - for (size_t regionIdx = 0; regionIdx < numRocktabTables; ++regionIdx) { - const RockwnodTable& rockwnodTable = rockwnodTables.template getTable(regionIdx); - const auto& rock2dtrTable = rock2dtrTables[regionIdx]; - - if (rockwnodTable.getSaturationColumn().size() != rock2dtrTable.sizeMultValues()) - throw std::runtime_error("Number of entries in ROCKWNOD and ROCK2DTR needs to match."); - - for (size_t xIdx = 0; xIdx < rock2dtrTable.size(); ++xIdx) { - rockCompTransMultWc_[regionIdx].appendXPos(rock2dtrTable.getPressureValue(xIdx)); - for (size_t yIdx = 0; yIdx < rockwnodTable.getSaturationColumn().size(); ++yIdx) - rockCompTransMultWc_[regionIdx].appendSamplePoint(xIdx, - rockwnodTable.getSaturationColumn()[yIdx], - rock2dtrTable.getTransMultValue(xIdx, yIdx)); - } - } - } - } -} - -template -Scalar EclGenericProblem:: -rockCompressibility(unsigned globalSpaceIdx) const -{ - if (this->rockParams_.empty()) - return 0.0; - - unsigned tableIdx = 0; - if (!this->rockTableIdx_.empty()) { - tableIdx = this->rockTableIdx_[globalSpaceIdx]; - } - return this->rockParams_[tableIdx].compressibility; -} - -template -Scalar EclGenericProblem:: -rockReferencePressure(unsigned globalSpaceIdx) const -{ - if (this->rockParams_.empty()) - return 1e5; - - unsigned tableIdx = 0; - if (!this->rockTableIdx_.empty()) { - tableIdx = this->rockTableIdx_[globalSpaceIdx]; - } - return this->rockParams_[tableIdx].referencePressure; -} - -template -Scalar EclGenericProblem:: -porosity(unsigned globalSpaceIdx, unsigned timeIdx) const -{ - return this->referencePorosity_[timeIdx][globalSpaceIdx]; -} - -template -Scalar EclGenericProblem:: -rockFraction(unsigned elementIdx, unsigned timeIdx) const -{ - const auto& fp = eclState_.fieldProps(); - const std::vector& poro = fp.get_double("PORO"); - // the reference porosity is defined as the accumulated pore volume divided by the - // geometric volume of the element. Note that it can - // be larger than 1.0 if porevolume multipliers are used - // to for instance implement larger boundary cells - Scalar porosity = poro[elementIdx]; - return referencePorosity(elementIdx, timeIdx) / porosity * (1 - porosity); -} - -template -template -void EclGenericProblem:: -updateNum(const std::string& name, std::vector& numbers, size_t num_regions) -{ - if (!eclState_.fieldProps().has_int(name)) - return; - - const auto& numData = eclState_.fieldProps().get_int(name); - - unsigned numElems = gridView_.size(/*codim=*/0); - numbers.resize(numElems); - for (unsigned elemIdx = 0; elemIdx < numElems; ++elemIdx) { - if (numData[elemIdx] > (int)num_regions) { - throw std::runtime_error("Values larger than maximum number of regions " + std::to_string(num_regions) + " provided in " + name); - } else if (numData[elemIdx] > 0) { - numbers[elemIdx] = static_cast(numData[elemIdx]) - 1; - } else { - throw std::runtime_error("zero or negative values provided for region array: " + name); - } - } -} - -template -void EclGenericProblem:: -updatePvtnum_() -{ - const auto num_regions = eclState_.getTableManager().getTabdims().getNumPVTTables(); - updateNum("PVTNUM", pvtnum_, num_regions); -} - -template -void EclGenericProblem:: -updateSatnum_() -{ - const auto num_regions = eclState_.getTableManager().getTabdims().getNumSatTables(); - updateNum("SATNUM", satnum_, num_regions); -} - -template -void EclGenericProblem:: -updateMiscnum_() -{ - const auto num_regions = 1; // we only support single region - updateNum("MISCNUM", miscnum_, num_regions); -} - -template -void EclGenericProblem:: -updatePlmixnum_() -{ - const auto num_regions = 1; // we only support single region - updateNum("PLMIXNUM", plmixnum_, num_regions); -} - -template -void EclGenericProblem:: -updateKrnum_() -{ - const auto num_regions = eclState_.getTableManager().getTabdims().getNumSatTables(); - updateNum("KRNUMX", krnumx_, num_regions); - updateNum("KRNUMY", krnumy_, num_regions); - updateNum("KRNUMZ", krnumz_, num_regions); -} - -template -void EclGenericProblem:: -updateImbnum_() -{ - const auto num_regions = eclState_.getTableManager().getTabdims().getNumSatTables(); - updateNum("IMBNUMX", imbnumx_, num_regions); - updateNum("IMBNUMY", imbnumy_, num_regions); - updateNum("IMBNUMZ", imbnumz_, num_regions); -} - -template -bool EclGenericProblem:: -vapparsActive(int episodeIdx) const -{ - const auto& oilVaporizationControl = schedule_[episodeIdx].oilvap(); - return (oilVaporizationControl.getType() == OilVaporizationProperties::OilVaporization::VAPPARS); -} - -template -bool EclGenericProblem:: -drsdtActive_(int episodeIdx) const -{ - const auto& oilVaporizationControl = schedule_[episodeIdx].oilvap(); - const bool bothOilGasActive = FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && - FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx); - return (oilVaporizationControl.drsdtActive() && bothOilGasActive); -} - -template -bool EclGenericProblem:: -drvdtActive_(int episodeIdx) const -{ - const auto& oilVaporizationControl = schedule_[episodeIdx].oilvap(); - const bool bothOilGasActive = FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && - FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx); - return (oilVaporizationControl.drvdtActive() && bothOilGasActive); -} - -template -bool EclGenericProblem:: -drsdtConvective_(int episodeIdx) const -{ - const auto& oilVaporizationControl = schedule_[episodeIdx].oilvap(); - const bool bothOilGasActive = FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && - FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx); - return (oilVaporizationControl.drsdtConvective() && bothOilGasActive); -} - -template -bool EclGenericProblem:: -beginEpisode_(bool enableExperiments, - int episodeIdx) -{ - if (enableExperiments && gridView_.comm().rank() == 0 && episodeIdx >= 0) { - // print some useful information in experimental mode. (the production - // simulator does this externally.) - std::ostringstream ss; - boost::posix_time::time_facet* facet = new boost::posix_time::time_facet("%d-%b-%Y"); - boost::posix_time::ptime curDateTime = - boost::posix_time::from_time_t(schedule_.simTime(episodeIdx)); - ss.imbue(std::locale(std::locale::classic(), facet)); - ss << "Report step " << episodeIdx + 1 - << "/" << schedule_.size() - 1 - << " at day " << schedule_.seconds(episodeIdx)/(24*3600) - << "/" << schedule_.seconds(schedule_.size() - 1)/(24*3600) - << ", date = " << curDateTime.date() - << "\n "; - OpmLog::info(ss.str()); - } - - const auto& events = schedule_[episodeIdx].events(); - - // react to TUNING changes - if (episodeIdx > 0 && enableTuning_ && events.hasEvent(ScheduleEvents::TUNING_CHANGE)) - { - const auto& sched_state = schedule_[episodeIdx]; - const auto& tuning = sched_state.tuning(); - initialTimeStepSize_ = sched_state.max_next_tstep(); - maxTimeStepAfterWellEvent_ = tuning.TMAXWC; - maxTimeStepSize_ = tuning.TSMAXZ; - restartShrinkFactor_ = 1./tuning.TSFCNV; - minTimeStepSize_ = tuning.TSMINZ; - return true; - } - - return false; -} - - -template -void EclGenericProblem:: -beginTimeStep_(bool enableExperiments, - int episodeIdx, - int timeStepIndex, - Scalar startTime, - Scalar time, - Scalar timeStepSize, - Scalar endTime) -{ - if (enableExperiments && gridView_.comm().rank() == 0 && episodeIdx >= 0) { - std::ostringstream ss; - boost::posix_time::time_facet* facet = new boost::posix_time::time_facet("%d-%b-%Y"); - boost::posix_time::ptime date = boost::posix_time::from_time_t(startTime) + boost::posix_time::milliseconds(static_cast(time / prefix::milli)); - ss.imbue(std::locale(std::locale::classic(), facet)); - ss <<"\nTime step " << timeStepIndex << ", stepsize " - << unit::convert::to(timeStepSize, unit::day) << " days," - << " at day " << (double)unit::convert::to(time, unit::day) - << "/" << (double)unit::convert::to(endTime, unit::day) - << ", date = " << date; - OpmLog::info(ss.str()); - } - - // update explicit quantities between timesteps. - const auto& oilVaporizationControl = schedule_[episodeIdx].oilvap(); - if (drsdtActive_(episodeIdx)) - // DRSDT is enabled - for (size_t pvtRegionIdx = 0; pvtRegionIdx < maxDRs_.size(); ++pvtRegionIdx) - maxDRs_[pvtRegionIdx] = oilVaporizationControl.getMaxDRSDT(pvtRegionIdx)*timeStepSize; - - if (drvdtActive_(episodeIdx)) - // DRVDT is enabled - for (size_t pvtRegionIdx = 0; pvtRegionIdx < maxDRv_.size(); ++pvtRegionIdx) - maxDRv_[pvtRegionIdx] = oilVaporizationControl.getMaxDRVDT(pvtRegionIdx)*timeStepSize; -} - -template -void EclGenericProblem:: -initFluidSystem_() -{ - FluidSystem::initFromState(eclState_, schedule_); -} - -template -void EclGenericProblem:: -readBlackoilExtentionsInitialConditions_(size_t numDof, - bool enableSolvent, - bool enablePolymer, - bool enablePolymerMolarWeight, - bool enableMICP) -{ - if (enableSolvent) { - if (eclState_.fieldProps().has_double("SSOL")) - solventSaturation_ = eclState_.fieldProps().get_double("SSOL"); - else - solventSaturation_.resize(numDof, 0.0); - } - - if (enablePolymer) { - if (eclState_.fieldProps().has_double("SPOLY")) - polymer_.concentration = eclState_.fieldProps().get_double("SPOLY"); - else - polymer_.concentration.resize(numDof, 0.0); - } - - if (enablePolymerMolarWeight) { - if (eclState_.fieldProps().has_double("SPOLYMW")) - polymer_.moleWeight = eclState_.fieldProps().get_double("SPOLYMW"); - else - polymer_.moleWeight.resize(numDof, 0.0); - } - - if (enableMICP) { - if (eclState_.fieldProps().has_double("SMICR")) - micp_.microbialConcentration = eclState_.fieldProps().get_double("SMICR"); - else - micp_.microbialConcentration.resize(numDof, 0.0); - if (eclState_.fieldProps().has_double("SOXYG")) - micp_.oxygenConcentration = eclState_.fieldProps().get_double("SOXYG"); - else - micp_.oxygenConcentration.resize(numDof, 0.0); - if (eclState_.fieldProps().has_double("SUREA")) - micp_.ureaConcentration = eclState_.fieldProps().get_double("SUREA"); - else - micp_.ureaConcentration.resize(numDof, 0.0); - if (eclState_.fieldProps().has_double("SBIOF")) - micp_.biofilmConcentration = eclState_.fieldProps().get_double("SBIOF"); - else - micp_.biofilmConcentration.resize(numDof, 0.0); - if (eclState_.fieldProps().has_double("SCALC")) - micp_.calciteConcentration = eclState_.fieldProps().get_double("SCALC"); - else - micp_.calciteConcentration.resize(numDof, 0.0); -} -} - - -template -Scalar EclGenericProblem:: -maxWaterSaturation(unsigned globalDofIdx) const -{ - if (maxWaterSaturation_.empty()) - return 0.0; - - return maxWaterSaturation_[globalDofIdx]; -} - -template -Scalar EclGenericProblem:: -minOilPressure(unsigned globalDofIdx) const -{ - if (minOilPressure_.empty()) - return 0.0; - - return minOilPressure_[globalDofIdx]; -} - -template -Scalar EclGenericProblem:: -overburdenPressure(unsigned elementIdx) const -{ - if (overburdenPressure_.empty()) - return 0.0; - - return overburdenPressure_[elementIdx]; -} - -template -Scalar EclGenericProblem:: -solventSaturation(unsigned elemIdx) const -{ - if (solventSaturation_.empty()) - return 0; - - return solventSaturation_[elemIdx]; -} - -template -Scalar EclGenericProblem:: -drsdtcon(unsigned elemIdx, int episodeIdx) const -{ - if (convectiveDrs_.empty()) - return 0; - - // The episode index is set to -1 in the initialization phase. - // Output drsdt value for index 0 - episodeIdx = std::max(episodeIdx, 0); - const auto& oilVaporizationControl = schedule_[episodeIdx].oilvap(); - int pvtRegionIdx = pvtRegionIndex(elemIdx); - return oilVaporizationControl.getMaxDRSDT(pvtRegionIdx)*convectiveDrs_[elemIdx]; -} - -template -Scalar EclGenericProblem:: -polymerConcentration(unsigned elemIdx) const -{ - if (polymer_.concentration.empty()) { - return 0; - } - - return polymer_.concentration[elemIdx]; -} - -template -Scalar EclGenericProblem:: -polymerMolecularWeight(const unsigned elemIdx) const -{ - if (polymer_.moleWeight.empty()) { - return 0.0; - } - - return polymer_.moleWeight[elemIdx]; -} - -template -Scalar EclGenericProblem:: -microbialConcentration(unsigned elemIdx) const -{ - if (micp_.microbialConcentration.empty()) { - return 0; - } - - return micp_.microbialConcentration[elemIdx]; -} - -template -Scalar EclGenericProblem:: -oxygenConcentration(unsigned elemIdx) const -{ - if (micp_.oxygenConcentration.empty()) { - return 0; - } - - return micp_.oxygenConcentration[elemIdx]; -} - -template -Scalar EclGenericProblem:: -ureaConcentration(unsigned elemIdx) const -{ - if (micp_.ureaConcentration.empty()) { - return 0; - } - - return micp_.ureaConcentration[elemIdx]; -} - -template -Scalar EclGenericProblem:: -biofilmConcentration(unsigned elemIdx) const -{ - if (micp_.biofilmConcentration.empty()) { - return 0; - } - - return micp_.biofilmConcentration[elemIdx]; -} - -template -Scalar EclGenericProblem:: -calciteConcentration(unsigned elemIdx) const -{ - if (micp_.calciteConcentration.empty()) { - return 0; - } - - return micp_.calciteConcentration[elemIdx]; -} - -template -unsigned EclGenericProblem:: -pvtRegionIndex(unsigned elemIdx) const -{ - if (pvtnum_.empty()) - return 0; - - return pvtnum_[elemIdx]; -} - -template -unsigned EclGenericProblem:: -satnumRegionIndex(unsigned elemIdx) const -{ - if (satnum_.empty()) - return 0; - - return satnum_[elemIdx]; -} - -template -unsigned EclGenericProblem:: -miscnumRegionIndex(unsigned elemIdx) const -{ - if (miscnum_.empty()) - return 0; - - return miscnum_[elemIdx]; -} - -template -unsigned EclGenericProblem:: -plmixnumRegionIndex(unsigned elemIdx) const -{ - if (plmixnum_.empty()) - return 0; - - return plmixnum_[elemIdx]; -} - -template -Scalar EclGenericProblem:: -maxPolymerAdsorption(unsigned elemIdx) const -{ - if (polymer_.maxAdsorption.empty()) { - return 0; - } - - return polymer_.maxAdsorption[elemIdx]; -} - -template -void EclGenericProblem:: -initDRSDT_(size_t numDof, - int episodeIdx) -{ - // deal with DRSDT - unsigned ntpvt = eclState_.runspec().tabdims().getNumPVTTables(); - //TODO We may want to only allocate these properties only if active. - //But since they may be activated at later time we need some more - //intrastructure to handle it - if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) { - maxDRv_.resize(ntpvt, 1e30); - lastRv_.resize(numDof, 0.0); - maxDRs_.resize(ntpvt, 1e30); - dRsDtOnlyFreeGas_.resize(ntpvt, false); - lastRs_.resize(numDof, 0.0); - maxDRv_.resize(ntpvt, 1e30); - lastRv_.resize(numDof, 0.0); - maxOilSaturation_.resize(numDof, 0.0); - if (drsdtConvective_(episodeIdx)) { - convectiveDrs_.resize(numDof, 1.0); - } - } -} - -template -bool EclGenericProblem:: -operator==(const EclGenericProblem& rhs) const -{ - return this->maxOilSaturation_ == rhs.maxOilSaturation_ && - this->maxWaterSaturation_ == rhs.maxWaterSaturation_ && - this->minOilPressure_ == rhs.minOilPressure_ && - this->overburdenPressure_ == rhs.overburdenPressure_ && - this->polymer_ == rhs.polymer_ && - this->solventSaturation_ == rhs.solventSaturation_ && - this->micp_ == rhs.micp_ && - this->lastRv_ == rhs.lastRv_ && - this->maxDRv_ == rhs.maxDRv_ && - this->convectiveDrs_ == rhs.convectiveDrs_ && - this->lastRs_ == rhs.lastRs_ && - this->maxDRs_ == rhs.maxDRs_ && - this->dRsDtOnlyFreeGas_ == rhs.dRsDtOnlyFreeGas_; -} - #if HAVE_DUNE_FEM template class EclGenericProblem>>, BlackOilFluidSystem, double>; template class EclGenericProblem>, - Opm::BlackOilFluidSystem< - double, - Opm::BlackOilDefaultIndexTraits>, - double>; -#if HAVE_DUNE_ALUGRID - -#if HAVE_MPI - using ALUGrid3CN = Dune::ALUGrid<3, 3, Dune::cube, Dune::nonconforming, Dune::ALUGridMPIComm>; -#else - using ALUGrid3CN = Dune::ALUGrid<3, 3, Dune::cube, Dune::nonconforming, Dune::ALUGridNoComm>; -#endif //HAVE_MPI - -template class EclGenericProblem>>, + Dune::Fem::AdaptiveLeafGridPart< + Dune::CpGrid, + Dune::PartitionIteratorType(4), + false> >, BlackOilFluidSystem, double>; -template class EclGenericProblem>, - Opm::BlackOilFluidSystem< - double, - Opm::BlackOilDefaultIndexTraits>, - double>; - -#endif //HAVE_DUNE_ALUGRID #else template class EclGenericProblem>, BlackOilFluidSystem, double>; -#if HAVE_DUNE_ALUGRID -#if HAVE_MPI - using ALUGrid3CN = const Dune::ALUGrid<3, 3, Dune::cube, Dune::nonconforming, Dune::ALUGridMPIComm>; -#else - using ALUGrid3CN = const Dune::ALUGrid<3, 3, Dune::cube, Dune::nonconforming, Dune::ALUGridNoComm>; -#endif //HAVE_MPI +#endif -template class EclGenericProblem>, - BlackOilFluidSystem, - double>; - -#endif //HAVE_DUNE_ALUGRID -#endif //HAVE_DUNE_FEM - -template class EclGenericProblem>, - BlackOilFluidSystem, - double>; - -} // namespace Opm +} // end namespace Opm diff --git a/ebos/eclgenericproblem_impl.hh b/ebos/eclgenericproblem_impl.hh new file mode 100644 index 000000000..479b434a1 --- /dev/null +++ b/ebos/eclgenericproblem_impl.hh @@ -0,0 +1,829 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/* + 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 2 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 . + + Consult the COPYING file in the top-level source directory of this + module for the precise wording of the license and the list of + copyright holders. +*/ +#ifndef EWOMS_GENERIC_ECL_PROBLEM_IMPL_HH +#define EWOMS_GENERIC_ECL_PROBLEM_IPML_HH + +#include + +#include + +#include + +#include +#include +#include +#include +#include + +#include + +#include +#include +#include +#include + +namespace Opm { + +int eclPositionalParameter(Dune::ParameterTree& tree, + std::set& seenParams, + std::string& errorMsg, + const char** argv, + int paramIdx) +{ + std::string param = argv[paramIdx]; + size_t i = param.find('='); + if (i != std::string::npos) { + std::string oldParamName = param.substr(0, i); + std::string oldParamValue = param.substr(i+1); + std::string newParamName = "--" + oldParamName; + std::replace(newParamName.begin(), + newParamName.end(), '_' , '-'); + errorMsg = + "The old syntax to specify parameters on the command line is no longer supported: " + "Try replacing '" + oldParamName + "=" + oldParamValue + "' with "+ + "'" + newParamName + "=" + oldParamValue + "'!"; + return 0; + } + + if (seenParams.count("EclDeckFileName") > 0) { + errorMsg = + "Parameter 'EclDeckFileName' specified multiple times" + " as a command line parameter"; + return 0; + } + + tree["EclDeckFileName"] = argv[paramIdx]; + seenParams.insert("EclDeckFileName"); + return 1; +} + +template +EclGenericProblem:: +EclGenericProblem(const EclipseState& eclState, + const Schedule& schedule, + const GridView& gridView) + : eclState_(eclState) + , schedule_(schedule) + , gridView_(gridView) +{ +} + +template +EclGenericProblem +EclGenericProblem:: +serializationTestObject(const EclipseState& eclState, + const Schedule& schedule, + const GridView& gridView) +{ + EclGenericProblem result(eclState, schedule, gridView); + result.maxOilSaturation_ = {1.0, 2.0}; + result.maxWaterSaturation_ = {6.0}; + result.minOilPressure_ = {7.0, 8.0, 9.0, 10.0}; + result.overburdenPressure_ = {11.0}; + result.solventSaturation_ = {15.0}; + result.polymer_ = PolymerSolutionContainer::serializationTestObject(); + result.micp_ = MICPSolutionContainer::serializationTestObject(); + result.lastRv_ = {21.0}; + result.maxDRv_ = {22.0, 23.0}; + result.convectiveDrs_ = {24.0, 25.0, 26.0}; + result.lastRs_ = {27.0}; + result.maxDRs_ = {28.0}; + result.dRsDtOnlyFreeGas_ = {false, true}; + + return result; +} + +template +std::string +EclGenericProblem:: +helpPreamble(int, + const char **argv) +{ + std::string desc = EclGenericProblem::briefDescription(); + if (!desc.empty()) + desc = desc + "\n"; + + return + "Usage: "+std::string(argv[0]) + " [OPTIONS] [ECL_DECK_FILENAME]\n" + + desc; +} + +template +std::string +EclGenericProblem:: +briefDescription() +{ + if (briefDescription_.empty()) + return + "The Ecl-deck Black-Oil reservoir Simulator (ebos); a hydrocarbon " + "reservoir simulation program that processes ECL-formatted input " + "files that is part of the Open Porous Media project " + "(https://opm-project.org).\n" + "\n" + "THE GOAL OF THE `ebos` SIMULATOR IS TO CATER FOR THE NEEDS OF " + "DEVELOPMENT AND RESEARCH. No guarantees are made for production use!"; + else + return briefDescription_; +} + +template +void EclGenericProblem:: +readRockParameters_(const std::vector& cellCenterDepths) +{ + const auto& rock_config = eclState_.getSimulationConfig().rock_config(); + + // read the rock compressibility parameters + { + const auto& comp = rock_config.comp(); + rockParams_.clear(); + for (const auto& c : comp) + rockParams_.push_back( { c.pref, c.compressibility } ); + } + + // read the parameters for water-induced rock compaction + readRockCompactionParameters_(); + + unsigned numElem = gridView_.size(0); + if (eclState_.fieldProps().has_int(rock_config.rocknum_property())) { + rockTableIdx_.resize(numElem); + const auto& num = eclState_.fieldProps().get_int(rock_config.rocknum_property()); + for (size_t elemIdx = 0; elemIdx < numElem; ++ elemIdx) { + rockTableIdx_[elemIdx] = num[elemIdx] - 1; + } + } + + // Store overburden pressure pr element + const auto& overburdTables = eclState_.getTableManager().getOverburdTables(); + if (!overburdTables.empty()) { + overburdenPressure_.resize(numElem,0.0); + size_t numRocktabTables = rock_config.num_rock_tables(); + + if (overburdTables.size() != numRocktabTables) + throw std::runtime_error(std::to_string(numRocktabTables) +" OVERBURD tables is expected, but " + std::to_string(overburdTables.size()) +" is provided"); + + std::vector> overburdenTables(numRocktabTables); + for (size_t regionIdx = 0; regionIdx < numRocktabTables; ++regionIdx) { + const OverburdTable& overburdTable = overburdTables.template getTable(regionIdx); + overburdenTables[regionIdx].setXYContainers(overburdTable.getDepthColumn(),overburdTable.getOverburdenPressureColumn()); + } + + for (size_t elemIdx = 0; elemIdx < numElem; ++ elemIdx) { + unsigned tableIdx = 0; + if (!rockTableIdx_.empty()) { + tableIdx = rockTableIdx_[elemIdx]; + } + overburdenPressure_[elemIdx] = overburdenTables[tableIdx].eval(cellCenterDepths[elemIdx], /*extrapolation=*/true); + } + } +} + +template +void EclGenericProblem:: +readRockCompactionParameters_() +{ + const auto& rock_config = eclState_.getSimulationConfig().rock_config(); + + if (!rock_config.active()) + return; // deck does not enable rock compaction + + unsigned numElem = gridView_.size(0); + switch (rock_config.hysteresis_mode()) { + case RockConfig::Hysteresis::REVERS: + break; + case RockConfig::Hysteresis::IRREVERS: + // interpolate the porv volume multiplier using the minimum pressure in the cell + // i.e. don't allow re-inflation. + minOilPressure_.resize(numElem, 1e99); + break; + default: + throw std::runtime_error("Not support ROCKOMP hysteresis option "); + } + + size_t numRocktabTables = rock_config.num_rock_tables(); + bool waterCompaction = rock_config.water_compaction(); + + if (!waterCompaction) { + const auto& rocktabTables = eclState_.getTableManager().getRocktabTables(); + if (rocktabTables.size() != numRocktabTables) + throw std::runtime_error("ROCKCOMP is activated." + std::to_string(numRocktabTables) + +" ROCKTAB tables is expected, but " + std::to_string(rocktabTables.size()) +" is provided"); + + rockCompPoroMult_.resize(numRocktabTables); + rockCompTransMult_.resize(numRocktabTables); + for (size_t regionIdx = 0; regionIdx < numRocktabTables; ++regionIdx) { + const auto& rocktabTable = rocktabTables.template getTable(regionIdx); + const auto& pressureColumn = rocktabTable.getPressureColumn(); + const auto& poroColumn = rocktabTable.getPoreVolumeMultiplierColumn(); + const auto& transColumn = rocktabTable.getTransmissibilityMultiplierColumn(); + rockCompPoroMult_[regionIdx].setXYContainers(pressureColumn, poroColumn); + rockCompTransMult_[regionIdx].setXYContainers(pressureColumn, transColumn); + } + } else { + const auto& rock2dTables = eclState_.getTableManager().getRock2dTables(); + const auto& rock2dtrTables = eclState_.getTableManager().getRock2dtrTables(); + const auto& rockwnodTables = eclState_.getTableManager().getRockwnodTables(); + maxWaterSaturation_.resize(numElem, 0.0); + + if (rock2dTables.size() != numRocktabTables) + throw std::runtime_error("Water compation option is selected in ROCKCOMP." + std::to_string(numRocktabTables) + +" ROCK2D tables is expected, but " + std::to_string(rock2dTables.size()) +" is provided"); + + if (rockwnodTables.size() != numRocktabTables) + throw std::runtime_error("Water compation option is selected in ROCKCOMP." + std::to_string(numRocktabTables) + +" ROCKWNOD tables is expected, but " + std::to_string(rockwnodTables.size()) +" is provided"); + //TODO check size match + rockCompPoroMultWc_.resize(numRocktabTables, TabulatedTwoDFunction(TabulatedTwoDFunction::InterpolationPolicy::Vertical)); + for (size_t regionIdx = 0; regionIdx < numRocktabTables; ++regionIdx) { + const RockwnodTable& rockwnodTable = rockwnodTables.template getTable(regionIdx); + const auto& rock2dTable = rock2dTables[regionIdx]; + + if (rockwnodTable.getSaturationColumn().size() != rock2dTable.sizeMultValues()) + throw std::runtime_error("Number of entries in ROCKWNOD and ROCK2D needs to match."); + + for (size_t xIdx = 0; xIdx < rock2dTable.size(); ++xIdx) { + rockCompPoroMultWc_[regionIdx].appendXPos(rock2dTable.getPressureValue(xIdx)); + for (size_t yIdx = 0; yIdx < rockwnodTable.getSaturationColumn().size(); ++yIdx) + rockCompPoroMultWc_[regionIdx].appendSamplePoint(xIdx, + rockwnodTable.getSaturationColumn()[yIdx], + rock2dTable.getPvmultValue(xIdx, yIdx)); + } + } + + if (!rock2dtrTables.empty()) { + rockCompTransMultWc_.resize(numRocktabTables, TabulatedTwoDFunction(TabulatedTwoDFunction::InterpolationPolicy::Vertical)); + for (size_t regionIdx = 0; regionIdx < numRocktabTables; ++regionIdx) { + const RockwnodTable& rockwnodTable = rockwnodTables.template getTable(regionIdx); + const auto& rock2dtrTable = rock2dtrTables[regionIdx]; + + if (rockwnodTable.getSaturationColumn().size() != rock2dtrTable.sizeMultValues()) + throw std::runtime_error("Number of entries in ROCKWNOD and ROCK2DTR needs to match."); + + for (size_t xIdx = 0; xIdx < rock2dtrTable.size(); ++xIdx) { + rockCompTransMultWc_[regionIdx].appendXPos(rock2dtrTable.getPressureValue(xIdx)); + for (size_t yIdx = 0; yIdx < rockwnodTable.getSaturationColumn().size(); ++yIdx) + rockCompTransMultWc_[regionIdx].appendSamplePoint(xIdx, + rockwnodTable.getSaturationColumn()[yIdx], + rock2dtrTable.getTransMultValue(xIdx, yIdx)); + } + } + } + } +} + +template +Scalar EclGenericProblem:: +rockCompressibility(unsigned globalSpaceIdx) const +{ + if (this->rockParams_.empty()) + return 0.0; + + unsigned tableIdx = 0; + if (!this->rockTableIdx_.empty()) { + tableIdx = this->rockTableIdx_[globalSpaceIdx]; + } + return this->rockParams_[tableIdx].compressibility; +} + +template +Scalar EclGenericProblem:: +rockReferencePressure(unsigned globalSpaceIdx) const +{ + if (this->rockParams_.empty()) + return 1e5; + + unsigned tableIdx = 0; + if (!this->rockTableIdx_.empty()) { + tableIdx = this->rockTableIdx_[globalSpaceIdx]; + } + return this->rockParams_[tableIdx].referencePressure; +} + +template +Scalar EclGenericProblem:: +porosity(unsigned globalSpaceIdx, unsigned timeIdx) const +{ + return this->referencePorosity_[timeIdx][globalSpaceIdx]; +} + +template +Scalar EclGenericProblem:: +rockFraction(unsigned elementIdx, unsigned timeIdx) const +{ + const auto& fp = eclState_.fieldProps(); + const std::vector& poro = fp.get_double("PORO"); + // the reference porosity is defined as the accumulated pore volume divided by the + // geometric volume of the element. Note that it can + // be larger than 1.0 if porevolume multipliers are used + // to for instance implement larger boundary cells + Scalar porosity = poro[elementIdx]; + return referencePorosity(elementIdx, timeIdx) / porosity * (1 - porosity); +} + +template +template +void EclGenericProblem:: +updateNum(const std::string& name, std::vector& numbers, size_t num_regions) +{ + if (!eclState_.fieldProps().has_int(name)) + return; + + const auto& numData = eclState_.fieldProps().get_int(name); + + unsigned numElems = gridView_.size(/*codim=*/0); + numbers.resize(numElems); + for (unsigned elemIdx = 0; elemIdx < numElems; ++elemIdx) { + if (numData[elemIdx] > (int)num_regions) { + throw std::runtime_error("Values larger than maximum number of regions " + std::to_string(num_regions) + " provided in " + name); + } else if (numData[elemIdx] > 0) { + numbers[elemIdx] = static_cast(numData[elemIdx]) - 1; + } else { + throw std::runtime_error("zero or negative values provided for region array: " + name); + } + } +} + +template +void EclGenericProblem:: +updatePvtnum_() +{ + const auto num_regions = eclState_.getTableManager().getTabdims().getNumPVTTables(); + updateNum("PVTNUM", pvtnum_, num_regions); +} + +template +void EclGenericProblem:: +updateSatnum_() +{ + const auto num_regions = eclState_.getTableManager().getTabdims().getNumSatTables(); + updateNum("SATNUM", satnum_, num_regions); +} + +template +void EclGenericProblem:: +updateMiscnum_() +{ + const auto num_regions = 1; // we only support single region + updateNum("MISCNUM", miscnum_, num_regions); +} + +template +void EclGenericProblem:: +updatePlmixnum_() +{ + const auto num_regions = 1; // we only support single region + updateNum("PLMIXNUM", plmixnum_, num_regions); +} + +template +void EclGenericProblem:: +updateKrnum_() +{ + const auto num_regions = eclState_.getTableManager().getTabdims().getNumSatTables(); + updateNum("KRNUMX", krnumx_, num_regions); + updateNum("KRNUMY", krnumy_, num_regions); + updateNum("KRNUMZ", krnumz_, num_regions); +} + +template +void EclGenericProblem:: +updateImbnum_() +{ + const auto num_regions = eclState_.getTableManager().getTabdims().getNumSatTables(); + updateNum("IMBNUMX", imbnumx_, num_regions); + updateNum("IMBNUMY", imbnumy_, num_regions); + updateNum("IMBNUMZ", imbnumz_, num_regions); +} + +template +bool EclGenericProblem:: +vapparsActive(int episodeIdx) const +{ + const auto& oilVaporizationControl = schedule_[episodeIdx].oilvap(); + return (oilVaporizationControl.getType() == OilVaporizationProperties::OilVaporization::VAPPARS); +} + +template +bool EclGenericProblem:: +drsdtActive_(int episodeIdx) const +{ + const auto& oilVaporizationControl = schedule_[episodeIdx].oilvap(); + const bool bothOilGasActive = FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && + FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx); + return (oilVaporizationControl.drsdtActive() && bothOilGasActive); +} + +template +bool EclGenericProblem:: +drvdtActive_(int episodeIdx) const +{ + const auto& oilVaporizationControl = schedule_[episodeIdx].oilvap(); + const bool bothOilGasActive = FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && + FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx); + return (oilVaporizationControl.drvdtActive() && bothOilGasActive); +} + +template +bool EclGenericProblem:: +drsdtConvective_(int episodeIdx) const +{ + const auto& oilVaporizationControl = schedule_[episodeIdx].oilvap(); + const bool bothOilGasActive = FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && + FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx); + return (oilVaporizationControl.drsdtConvective() && bothOilGasActive); +} + +template +bool EclGenericProblem:: +beginEpisode_(bool enableExperiments, + int episodeIdx) +{ + if (enableExperiments && gridView_.comm().rank() == 0 && episodeIdx >= 0) { + // print some useful information in experimental mode. (the production + // simulator does this externally.) + std::ostringstream ss; + boost::posix_time::time_facet* facet = new boost::posix_time::time_facet("%d-%b-%Y"); + boost::posix_time::ptime curDateTime = + boost::posix_time::from_time_t(schedule_.simTime(episodeIdx)); + ss.imbue(std::locale(std::locale::classic(), facet)); + ss << "Report step " << episodeIdx + 1 + << "/" << schedule_.size() - 1 + << " at day " << schedule_.seconds(episodeIdx)/(24*3600) + << "/" << schedule_.seconds(schedule_.size() - 1)/(24*3600) + << ", date = " << curDateTime.date() + << "\n "; + OpmLog::info(ss.str()); + } + + const auto& events = schedule_[episodeIdx].events(); + + // react to TUNING changes + if (episodeIdx > 0 && enableTuning_ && events.hasEvent(ScheduleEvents::TUNING_CHANGE)) + { + const auto& sched_state = schedule_[episodeIdx]; + const auto& tuning = sched_state.tuning(); + initialTimeStepSize_ = sched_state.max_next_tstep(); + maxTimeStepAfterWellEvent_ = tuning.TMAXWC; + maxTimeStepSize_ = tuning.TSMAXZ; + restartShrinkFactor_ = 1./tuning.TSFCNV; + minTimeStepSize_ = tuning.TSMINZ; + return true; + } + + return false; +} + + +template +void EclGenericProblem:: +beginTimeStep_(bool enableExperiments, + int episodeIdx, + int timeStepIndex, + Scalar startTime, + Scalar time, + Scalar timeStepSize, + Scalar endTime) +{ + if (enableExperiments && gridView_.comm().rank() == 0 && episodeIdx >= 0) { + std::ostringstream ss; + boost::posix_time::time_facet* facet = new boost::posix_time::time_facet("%d-%b-%Y"); + boost::posix_time::ptime date = boost::posix_time::from_time_t(startTime) + boost::posix_time::milliseconds(static_cast(time / prefix::milli)); + ss.imbue(std::locale(std::locale::classic(), facet)); + ss <<"\nTime step " << timeStepIndex << ", stepsize " + << unit::convert::to(timeStepSize, unit::day) << " days," + << " at day " << (double)unit::convert::to(time, unit::day) + << "/" << (double)unit::convert::to(endTime, unit::day) + << ", date = " << date; + OpmLog::info(ss.str()); + } + + // update explicit quantities between timesteps. + const auto& oilVaporizationControl = schedule_[episodeIdx].oilvap(); + if (drsdtActive_(episodeIdx)) + // DRSDT is enabled + for (size_t pvtRegionIdx = 0; pvtRegionIdx < maxDRs_.size(); ++pvtRegionIdx) + maxDRs_[pvtRegionIdx] = oilVaporizationControl.getMaxDRSDT(pvtRegionIdx)*timeStepSize; + + if (drvdtActive_(episodeIdx)) + // DRVDT is enabled + for (size_t pvtRegionIdx = 0; pvtRegionIdx < maxDRv_.size(); ++pvtRegionIdx) + maxDRv_[pvtRegionIdx] = oilVaporizationControl.getMaxDRVDT(pvtRegionIdx)*timeStepSize; +} + +template +void EclGenericProblem:: +initFluidSystem_() +{ + FluidSystem::initFromState(eclState_, schedule_); +} + +template +void EclGenericProblem:: +readBlackoilExtentionsInitialConditions_(size_t numDof, + bool enableSolvent, + bool enablePolymer, + bool enablePolymerMolarWeight, + bool enableMICP) +{ + if (enableSolvent) { + if (eclState_.fieldProps().has_double("SSOL")) + solventSaturation_ = eclState_.fieldProps().get_double("SSOL"); + else + solventSaturation_.resize(numDof, 0.0); + } + + if (enablePolymer) { + if (eclState_.fieldProps().has_double("SPOLY")) { + polymer_.concentration = eclState_.fieldProps().get_double("SPOLY"); + } else { + polymer_.concentration.resize(numDof, 0.0); + } + } + + if (enablePolymerMolarWeight) { + if (eclState_.fieldProps().has_double("SPOLYMW")) { + polymer_.moleWeight = eclState_.fieldProps().get_double("SPOLYMW"); + } else { + polymer_.moleWeight.resize(numDof, 0.0); + } + } + + if (enableMICP) { + if (eclState_.fieldProps().has_double("SMICR")) { + micp_.microbialConcentration = eclState_.fieldProps().get_double("SMICR"); + } else { + micp_.microbialConcentration.resize(numDof, 0.0); + } + if (eclState_.fieldProps().has_double("SOXYG")) { + micp_.oxygenConcentration = eclState_.fieldProps().get_double("SOXYG"); + } else { + micp_.oxygenConcentration.resize(numDof, 0.0); + } + if (eclState_.fieldProps().has_double("SUREA")) { + micp_.ureaConcentration = eclState_.fieldProps().get_double("SUREA"); + } else { + micp_.ureaConcentration.resize(numDof, 0.0); + } + if (eclState_.fieldProps().has_double("SBIOF")) { + micp_.biofilmConcentration = eclState_.fieldProps().get_double("SBIOF"); + } else { + micp_.biofilmConcentration.resize(numDof, 0.0); + } + if (eclState_.fieldProps().has_double("SCALC")) { + micp_.calciteConcentration = eclState_.fieldProps().get_double("SCALC"); + } else { + micp_.calciteConcentration.resize(numDof, 0.0); + } +} +} + + +template +Scalar EclGenericProblem:: +maxWaterSaturation(unsigned globalDofIdx) const +{ + if (maxWaterSaturation_.empty()) + return 0.0; + + return maxWaterSaturation_[globalDofIdx]; +} + +template +Scalar EclGenericProblem:: +minOilPressure(unsigned globalDofIdx) const +{ + if (minOilPressure_.empty()) + return 0.0; + + return minOilPressure_[globalDofIdx]; +} + +template +Scalar EclGenericProblem:: +overburdenPressure(unsigned elementIdx) const +{ + if (overburdenPressure_.empty()) + return 0.0; + + return overburdenPressure_[elementIdx]; +} + +template +Scalar EclGenericProblem:: +solventSaturation(unsigned elemIdx) const +{ + if (solventSaturation_.empty()) + return 0; + + return solventSaturation_[elemIdx]; +} + +template +Scalar EclGenericProblem:: +drsdtcon(unsigned elemIdx, int episodeIdx) const +{ + if (convectiveDrs_.empty()) + return 0; + + // The episode index is set to -1 in the initialization phase. + // Output drsdt value for index 0 + episodeIdx = std::max(episodeIdx, 0); + const auto& oilVaporizationControl = schedule_[episodeIdx].oilvap(); + int pvtRegionIdx = pvtRegionIndex(elemIdx); + return oilVaporizationControl.getMaxDRSDT(pvtRegionIdx)*convectiveDrs_[elemIdx]; +} + +template +Scalar EclGenericProblem:: +polymerConcentration(unsigned elemIdx) const +{ + if (polymer_.concentration.empty()) { + return 0; + } + + return polymer_.concentration[elemIdx]; +} + +template +Scalar EclGenericProblem:: +polymerMolecularWeight(const unsigned elemIdx) const +{ + if (polymer_.moleWeight.empty()) { + return 0.0; + } + + return polymer_.moleWeight[elemIdx]; +} + +template +Scalar EclGenericProblem:: +microbialConcentration(unsigned elemIdx) const +{ + if (micp_.microbialConcentration.empty()) { + return 0; + } + + return micp_.microbialConcentration[elemIdx]; +} + +template +Scalar EclGenericProblem:: +oxygenConcentration(unsigned elemIdx) const +{ + if (micp_.oxygenConcentration.empty()) { + return 0; + } + + return micp_.oxygenConcentration[elemIdx]; +} + +template +Scalar EclGenericProblem:: +ureaConcentration(unsigned elemIdx) const +{ + if (micp_.ureaConcentration.empty()) { + return 0; + } + + return micp_.ureaConcentration[elemIdx]; +} + +template +Scalar EclGenericProblem:: +biofilmConcentration(unsigned elemIdx) const +{ + if (micp_.biofilmConcentration.empty()) { + return 0; + } + + return micp_.biofilmConcentration[elemIdx]; +} + +template +Scalar EclGenericProblem:: +calciteConcentration(unsigned elemIdx) const +{ + if (micp_.calciteConcentration.empty()) { + return 0; + } + + return micp_.calciteConcentration[elemIdx]; +} + +template +unsigned EclGenericProblem:: +pvtRegionIndex(unsigned elemIdx) const +{ + if (pvtnum_.empty()) + return 0; + + return pvtnum_[elemIdx]; +} + +template +unsigned EclGenericProblem:: +satnumRegionIndex(unsigned elemIdx) const +{ + if (satnum_.empty()) + return 0; + + return satnum_[elemIdx]; +} + +template +unsigned EclGenericProblem:: +miscnumRegionIndex(unsigned elemIdx) const +{ + if (miscnum_.empty()) + return 0; + + return miscnum_[elemIdx]; +} + +template +unsigned EclGenericProblem:: +plmixnumRegionIndex(unsigned elemIdx) const +{ + if (plmixnum_.empty()) + return 0; + + return plmixnum_[elemIdx]; +} + +template +Scalar EclGenericProblem:: +maxPolymerAdsorption(unsigned elemIdx) const +{ + if (polymer_.maxAdsorption.empty()) { + return 0; + } + + return polymer_.maxAdsorption[elemIdx]; +} + +template +void EclGenericProblem:: +initDRSDT_(size_t numDof, + int episodeIdx) +{ + // deal with DRSDT + unsigned ntpvt = eclState_.runspec().tabdims().getNumPVTTables(); + //TODO We may want to only allocate these properties only if active. + //But since they may be activated at later time we need some more + //intrastructure to handle it + if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) { + maxDRv_.resize(ntpvt, 1e30); + lastRv_.resize(numDof, 0.0); + maxDRs_.resize(ntpvt, 1e30); + dRsDtOnlyFreeGas_.resize(ntpvt, false); + lastRs_.resize(numDof, 0.0); + maxDRv_.resize(ntpvt, 1e30); + lastRv_.resize(numDof, 0.0); + maxOilSaturation_.resize(numDof, 0.0); + if (drsdtConvective_(episodeIdx)) { + convectiveDrs_.resize(numDof, 1.0); + } + } +} + +template +bool EclGenericProblem:: +operator==(const EclGenericProblem& rhs) const +{ + return this->maxOilSaturation_ == rhs.maxOilSaturation_ && + this->maxWaterSaturation_ == rhs.maxWaterSaturation_ && + this->minOilPressure_ == rhs.minOilPressure_ && + this->overburdenPressure_ == rhs.overburdenPressure_ && + this->solventSaturation_ == rhs.solventSaturation_ && + this->polymer_ == rhs.polymer_ && + this->micp_ == rhs.micp_ && + this->lastRv_ == rhs.lastRv_ && + this->maxDRv_ == rhs.maxDRv_ && + this->convectiveDrs_ == rhs.convectiveDrs_ && + this->lastRs_ == rhs.lastRs_ && + this->maxDRs_ == rhs.maxDRs_ && + this->dRsDtOnlyFreeGas_ == rhs.dRsDtOnlyFreeGas_; +} + +} // namespace Opm + +#endif // EWOMS_GENERIC_ECL_PROBLEM_IMPL_HH diff --git a/ebos/eclgenericthresholdpressure.cc b/ebos/eclgenericthresholdpressure.cc index 3d8c90547..1b18d252f 100644 --- a/ebos/eclgenericthresholdpressure.cc +++ b/ebos/eclgenericthresholdpressure.cc @@ -24,22 +24,9 @@ #include #include #include - -#include -#include -#include -#include -#include - -#include +#include #include -#include -#if HAVE_DUNE_ALUGRID -#include -#include -#include "alucartesianindexmapper.hh" -#endif // HAVE_DUNE_ALUGRID #if HAVE_DUNE_FEM #include @@ -47,272 +34,8 @@ #include #endif // HAVE_DUNE_FEM -#include - -#include -#include -#include - namespace Opm { -template -EclGenericThresholdPressure:: -EclGenericThresholdPressure(const CartesianIndexMapper& cartMapper, - const GridView& gridView, - const ElementMapper& elementMapper, - const EclipseState& eclState) - : cartMapper_(cartMapper) - , gridView_(gridView) - , elementMapper_(elementMapper) - , eclState_(eclState) -{ -} - -template -Scalar EclGenericThresholdPressure:: -thresholdPressure(int elem1Idx, int elem2Idx) const -{ - if (!enableThresholdPressure_) - return 0.0; - - // threshold pressure accross faults - if (!thpresftValues_.empty()) { - int cartElem1Idx = cartMapper_.cartesianIndex(elem1Idx); - int cartElem2Idx = cartMapper_.cartesianIndex(elem2Idx); - - assert(0 <= cartElem1Idx && static_cast(cartElemFaultIdx_.size()) > cartElem1Idx); - assert(0 <= cartElem2Idx && static_cast(cartElemFaultIdx_.size()) > cartElem2Idx); - - int fault1Idx = cartElemFaultIdx_[cartElem1Idx]; - int fault2Idx = cartElemFaultIdx_[cartElem2Idx]; - if (fault1Idx != -1 && fault1Idx == fault2Idx) - // inside a fault there's no threshold pressure, even accross EQUIL - // regions. - return 0.0; - if (fault1Idx != fault2Idx) { - // TODO: which value if a cell is part of multiple faults? we take - // the maximum here. - Scalar val1 = (fault1Idx >= 0) ? thpresftValues_[fault1Idx] : 0.0; - Scalar val2 = (fault2Idx >= 0) ? thpresftValues_[fault2Idx] : 0.0; - return std::max(val1, val2); - } - } - - // threshold pressure accross EQUIL regions - auto equilRegion1Idx = elemEquilRegion_[elem1Idx]; - auto equilRegion2Idx = elemEquilRegion_[elem2Idx]; - - if (equilRegion1Idx == equilRegion2Idx) - return 0.0; - - return thpres_[equilRegion1Idx*numEquilRegions_ + equilRegion2Idx]; -} - -template -void EclGenericThresholdPressure:: -finishInit() -{ - unsigned numElements = gridView_.size(/*codim=*/0); - - const auto& simConfig = eclState_.getSimulationConfig(); - - enableThresholdPressure_ = simConfig.useThresholdPressure(); - if (!enableThresholdPressure_) - return; - - numEquilRegions_ = eclState_.getTableManager().getEqldims().getNumEquilRegions(); - const decltype(numEquilRegions_) maxRegions = - std::numeric_limits>::max(); - - if (numEquilRegions_ > maxRegions) { - // make sure that the index of an equilibration region can be stored - // in the vector - OPM_THROW(std::invalid_argument, - (fmt::format("The maximum number of supported " - "equilibration regions by OPM flow is {}, but " - "{} are used!", - maxRegions, numEquilRegions_))); - } - - if (numEquilRegions_ > 2048) { - // warn about performance - OpmLog::warning(fmt::format("Number of equilibration regions is {}, which is " - "rather large. Note, that this might " - "have a negative impact on performance " - "and memory consumption.", numEquilRegions_)); - } - - // internalize the data specified using the EQLNUM keyword - const auto& fp = eclState_.fieldProps(); - const auto& equilRegionData = fp.get_int("EQLNUM"); - elemEquilRegion_.resize(numElements, 0); - for (unsigned elemIdx = 0; elemIdx < numElements; ++elemIdx) { - elemEquilRegion_[elemIdx] = equilRegionData[elemIdx] - 1; - } - - /* - If this is a restart run the ThresholdPressure object will be active, - but it will *not* be properly initialized with numerical values. The - values must instead come from the THPRES vector in the restart file. - */ - if (simConfig.getThresholdPressure().restart()) - return; - - // allocate the array which specifies the threshold pressures - thpres_.resize(numEquilRegions_*numEquilRegions_, 0.0); - thpresDefault_.resize(numEquilRegions_*numEquilRegions_, 0.0); -} - -template -void EclGenericThresholdPressure:: -applyExplicitThresholdPressures_() -{ - const SimulationConfig& simConfig = eclState_.getSimulationConfig(); - const auto& thpres = simConfig.getThresholdPressure(); - - // set the threshold pressures for all EQUIL region boundaries which have a - // intersection in the grid - for (const auto& elem : elements(gridView_, Dune::Partitions::interior)) { - for (const auto& intersection : intersections(gridView_, elem)) { - if (intersection.boundary()) - continue; // ignore boundary intersections for now (TODO?) - else if (!intersection.neighbor()) //processor boundary but not domain boundary - continue; - - const auto& inside = intersection.inside(); - const auto& outside = intersection.outside(); - - unsigned insideElemIdx = elementMapper_.index(inside); - unsigned outsideElemIdx = elementMapper_.index(outside); - - auto equilRegionInside = elemEquilRegion_[insideElemIdx]; - auto equilRegionOutside = elemEquilRegion_[outsideElemIdx]; - if (thpres.hasRegionBarrier(equilRegionInside + 1, equilRegionOutside + 1)) { - Scalar pth = 0.0; - if (thpres.hasThresholdPressure(equilRegionInside + 1, equilRegionOutside + 1)) { - // threshold pressure explicitly specified - pth = thpres.getThresholdPressure(equilRegionInside + 1, equilRegionOutside + 1); - } - else { - // take the threshold pressure from the initial condition - unsigned offset = equilRegionInside*numEquilRegions_ + equilRegionOutside; - pth = thpresDefault_[offset]; - } - - unsigned offset1 = equilRegionInside*numEquilRegions_ + equilRegionOutside; - unsigned offset2 = equilRegionOutside*numEquilRegions_ + equilRegionInside; - - thpres_[offset1] = pth; - thpres_[offset2] = pth; - } - } - } - - // apply threshold pressures across faults - if (thpres.ftSize() > 0) - configureThpresft_(); -} - -template -void EclGenericThresholdPressure:: -configureThpresft_() -{ - // retrieve the faults collection. - const FaultCollection& faults = eclState_.getFaults(); - const SimulationConfig& simConfig = eclState_.getSimulationConfig(); - const auto& thpres = simConfig.getThresholdPressure(); - - // extract the multipliers - int numFaults = faults.size(); - int numCartesianElem = eclState_.getInputGrid().getCartesianSize(); - thpresftValues_.resize(numFaults, -1.0); - cartElemFaultIdx_.resize(numCartesianElem, -1); - for (size_t faultIdx = 0; faultIdx < faults.size(); faultIdx++) { - auto& fault = faults.getFault(faultIdx); - thpresftValues_[faultIdx] = thpres.getThresholdPressureFault(faultIdx); - for (const FaultFace& face : fault) - // "face" is a misnomer because the object describes a set of cell - // indices, but we go with the conventions of the parser here... - for (size_t cartElemIdx : face) - cartElemFaultIdx_[cartElemIdx] = faultIdx; - } -} - -template -std::vector -EclGenericThresholdPressure:: -getRestartVector() const -{ - if (!enableThresholdPressure_) - return {}; - - std::vector result(numEquilRegions_ * numEquilRegions_, 0.0); - const auto& simConfig = eclState_.getSimulationConfig(); - const auto& thpres = simConfig.getThresholdPressure(); - - std::size_t idx = 0; - for (unsigned j = 1; j <= numEquilRegions_; ++j) { - for (unsigned i = 1; i <= numEquilRegions_; ++i, ++idx) { - if (thpres.hasRegionBarrier(i, j)) { - if (thpres.hasThresholdPressure(i, j)) { - result[idx] = thpres.getThresholdPressure(i, j); - } else { - result[idx] = this->thpresDefault_[idx]; - } - } - } - } - - return result; -} - -template -void -EclGenericThresholdPressure:: -logPressures() -{ - if (!enableThresholdPressure_) - return; - - auto lineFormat = [this](unsigned i, unsigned j, double val) - { - const auto& units = eclState_.getUnits(); - return fmt::format("{:4}{:>6}{:23}{:>6}{:24}{:>11.07}{:7}{}\n", - " ", i, - " ", j, - " ", units.from_si(UnitSystem::measure::pressure, val), - " ", units.name(UnitSystem::measure::pressure)); - }; - auto lineFormatS = [](auto s1, auto s2, auto s3) - { - return fmt::format("{:4}{:^16}{:13}{:^9}{:21}{:^18}\n", - " ", s1, " ", s2, " ", s3); - }; - - std::string str = "\nLIST OF ALL NON-ZERO THRESHOLD PRESSURES\n" - "----------------------------------------\n" - "\n"; - str += lineFormatS("FLOW FROM REGION", "TO REGION", "THRESHOLD PRESSURE"); - str += lineFormatS(std::string(16, '-'), std::string(9, '-'), std::string(18, '-')); - const auto& simConfig = eclState_.getSimulationConfig(); - const auto& thpres = simConfig.getThresholdPressure(); - - for (unsigned i = 1; i <= numEquilRegions_; ++i) { - for (unsigned j = (thpres.irreversible() ? 1 : i); j <= numEquilRegions_; ++j) { - if (thpres.hasRegionBarrier(i, j)) { - if (thpres.hasThresholdPressure(i, j)) { - str += lineFormat(i, j, thpres.getThresholdPressure(j, i)); - } else { - std::size_t idx = (j - 1) * numEquilRegions_ + (i - 1); - str += lineFormat(i, j, this->thpresDefault_[idx]); - } - } - } - } - str += lineFormatS(std::string(16, '-'), std::string(9, '-'), std::string(18, '-')); - OpmLog::note(str); -} - #if HAVE_DUNE_FEM template class EclGenericThresholdPressure>>, @@ -331,57 +54,12 @@ template class EclGenericThresholdPressure>>, double>; -#if HAVE_DUNE_ALUGRID -#if HAVE_MPI - using ALUGrid3CN = Dune::ALUGrid<3, 3, Dune::cube, Dune::nonconforming, Dune::ALUGridMPIComm>; -#else - using ALUGrid3CN = Dune::ALUGrid<3, 3, Dune::cube, Dune::nonconforming, Dune::ALUGridNoComm>; -#endif //HAVE_MPI - -template class EclGenericThresholdPressure>>, - Dune::MultipleCodimMultipleGeomTypeMapper>>>, - double>; -template class EclGenericThresholdPressure >, - Dune::MultipleCodimMultipleGeomTypeMapper< - Dune::Fem::GridPart2GridViewImpl< - Dune::Fem::AdaptiveLeafGridPart< - ALUGrid3CN, - Dune::PartitionIteratorType(4), - false>>>, - double>; - - - -#endif //HAVE_DUNE_ALUGRID #else template class EclGenericThresholdPressure>, Dune::MultipleCodimMultipleGeomTypeMapper>>, double>; -#if HAVE_DUNE_ALUGRID - -#if HAVE_MPI - using ALUGrid3CN = Dune::ALUGrid<3, 3, Dune::cube, Dune::nonconforming, Dune::ALUGridMPIComm>; -#else - using ALUGrid3CN = Dune::ALUGrid<3, 3, Dune::cube, Dune::nonconforming, Dune::ALUGridNoComm>; -#endif //HAVE_MPI -template class EclGenericThresholdPressure>, - Dune::MultipleCodimMultipleGeomTypeMapper>>, - double>; -#endif //HAVE_DUNE_ALUGRID - #endif -template class EclGenericThresholdPressure, - Dune::GridView>, - Dune::MultipleCodimMultipleGeomTypeMapper>>, - double>; } // namespace Opm diff --git a/ebos/eclgenericthresholdpressure_impl.hh b/ebos/eclgenericthresholdpressure_impl.hh new file mode 100644 index 000000000..bd9ac1d33 --- /dev/null +++ b/ebos/eclgenericthresholdpressure_impl.hh @@ -0,0 +1,307 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/* + 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 2 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 . + + Consult the COPYING file in the top-level source directory of this + module for the precise wording of the license and the list of + copyright holders. +*/ +#ifndef EWOMS_ECL_GENERIC_THRESHOLD_PRESSURE_IMPL_HH +#define EWOMS_ECL_GENERIC_THRESHOLD_PRESSURE_IMPL_HH + +#include + +#include +#include + +#include + +#include +#include +#include +#include +#include + +#include + +#include +#include +#include +#include + +namespace Opm { + +template +EclGenericThresholdPressure:: +EclGenericThresholdPressure(const CartesianIndexMapper& cartMapper, + const GridView& gridView, + const ElementMapper& elementMapper, + const EclipseState& eclState) + : cartMapper_(cartMapper) + , gridView_(gridView) + , elementMapper_(elementMapper) + , eclState_(eclState) +{ +} + +template +Scalar EclGenericThresholdPressure:: +thresholdPressure(int elem1Idx, int elem2Idx) const +{ + if (!enableThresholdPressure_) + return 0.0; + + // threshold pressure accross faults + if (!thpresftValues_.empty()) { + int cartElem1Idx = cartMapper_.cartesianIndex(elem1Idx); + int cartElem2Idx = cartMapper_.cartesianIndex(elem2Idx); + + assert(0 <= cartElem1Idx && static_cast(cartElemFaultIdx_.size()) > cartElem1Idx); + assert(0 <= cartElem2Idx && static_cast(cartElemFaultIdx_.size()) > cartElem2Idx); + + int fault1Idx = cartElemFaultIdx_[cartElem1Idx]; + int fault2Idx = cartElemFaultIdx_[cartElem2Idx]; + if (fault1Idx != -1 && fault1Idx == fault2Idx) + // inside a fault there's no threshold pressure, even accross EQUIL + // regions. + return 0.0; + if (fault1Idx != fault2Idx) { + // TODO: which value if a cell is part of multiple faults? we take + // the maximum here. + Scalar val1 = (fault1Idx >= 0) ? thpresftValues_[fault1Idx] : 0.0; + Scalar val2 = (fault2Idx >= 0) ? thpresftValues_[fault2Idx] : 0.0; + return std::max(val1, val2); + } + } + + // threshold pressure accross EQUIL regions + auto equilRegion1Idx = elemEquilRegion_[elem1Idx]; + auto equilRegion2Idx = elemEquilRegion_[elem2Idx]; + + if (equilRegion1Idx == equilRegion2Idx) + return 0.0; + + return thpres_[equilRegion1Idx*numEquilRegions_ + equilRegion2Idx]; +} + +template +void EclGenericThresholdPressure:: +finishInit() +{ + unsigned numElements = gridView_.size(/*codim=*/0); + + const auto& simConfig = eclState_.getSimulationConfig(); + + enableThresholdPressure_ = simConfig.useThresholdPressure(); + if (!enableThresholdPressure_) + return; + + numEquilRegions_ = eclState_.getTableManager().getEqldims().getNumEquilRegions(); + const decltype(numEquilRegions_) maxRegions = + std::numeric_limits>::max(); + + if (numEquilRegions_ > maxRegions) { + // make sure that the index of an equilibration region can be stored + // in the vector + OPM_THROW(std::invalid_argument, + (fmt::format("The maximum number of supported " + "equilibration regions by OPM flow is {}, but " + "{} are used!", + maxRegions, numEquilRegions_))); + } + + if (numEquilRegions_ > 2048) { + // warn about performance + OpmLog::warning(fmt::format("Number of equilibration regions is {}, which is " + "rather large. Note, that this might " + "have a negative impact on performance " + "and memory consumption.", numEquilRegions_)); + } + + // internalize the data specified using the EQLNUM keyword + const auto& fp = eclState_.fieldProps(); + const auto& equilRegionData = fp.get_int("EQLNUM"); + elemEquilRegion_.resize(numElements, 0); + for (unsigned elemIdx = 0; elemIdx < numElements; ++elemIdx) { + elemEquilRegion_[elemIdx] = equilRegionData[elemIdx] - 1; + } + + /* + If this is a restart run the ThresholdPressure object will be active, + but it will *not* be properly initialized with numerical values. The + values must instead come from the THPRES vector in the restart file. + */ + if (simConfig.getThresholdPressure().restart()) + return; + + // allocate the array which specifies the threshold pressures + thpres_.resize(numEquilRegions_*numEquilRegions_, 0.0); + thpresDefault_.resize(numEquilRegions_*numEquilRegions_, 0.0); +} + +template +void EclGenericThresholdPressure:: +applyExplicitThresholdPressures_() +{ + const SimulationConfig& simConfig = eclState_.getSimulationConfig(); + const auto& thpres = simConfig.getThresholdPressure(); + + // set the threshold pressures for all EQUIL region boundaries which have a + // intersection in the grid + for (const auto& elem : elements(gridView_, Dune::Partitions::interior)) { + for (const auto& intersection : intersections(gridView_, elem)) { + if (intersection.boundary()) + continue; // ignore boundary intersections for now (TODO?) + else if (!intersection.neighbor()) //processor boundary but not domain boundary + continue; + + const auto& inside = intersection.inside(); + const auto& outside = intersection.outside(); + + unsigned insideElemIdx = elementMapper_.index(inside); + unsigned outsideElemIdx = elementMapper_.index(outside); + + auto equilRegionInside = elemEquilRegion_[insideElemIdx]; + auto equilRegionOutside = elemEquilRegion_[outsideElemIdx]; + if (thpres.hasRegionBarrier(equilRegionInside + 1, equilRegionOutside + 1)) { + Scalar pth = 0.0; + if (thpres.hasThresholdPressure(equilRegionInside + 1, equilRegionOutside + 1)) { + // threshold pressure explicitly specified + pth = thpres.getThresholdPressure(equilRegionInside + 1, equilRegionOutside + 1); + } + else { + // take the threshold pressure from the initial condition + unsigned offset = equilRegionInside*numEquilRegions_ + equilRegionOutside; + pth = thpresDefault_[offset]; + } + + unsigned offset1 = equilRegionInside*numEquilRegions_ + equilRegionOutside; + unsigned offset2 = equilRegionOutside*numEquilRegions_ + equilRegionInside; + + thpres_[offset1] = pth; + thpres_[offset2] = pth; + } + } + } + + // apply threshold pressures across faults + if (thpres.ftSize() > 0) + configureThpresft_(); +} + +template +void EclGenericThresholdPressure:: +configureThpresft_() +{ + // retrieve the faults collection. + const FaultCollection& faults = eclState_.getFaults(); + const SimulationConfig& simConfig = eclState_.getSimulationConfig(); + const auto& thpres = simConfig.getThresholdPressure(); + + // extract the multipliers + int numFaults = faults.size(); + int numCartesianElem = eclState_.getInputGrid().getCartesianSize(); + thpresftValues_.resize(numFaults, -1.0); + cartElemFaultIdx_.resize(numCartesianElem, -1); + for (size_t faultIdx = 0; faultIdx < faults.size(); faultIdx++) { + auto& fault = faults.getFault(faultIdx); + thpresftValues_[faultIdx] = thpres.getThresholdPressureFault(faultIdx); + for (const FaultFace& face : fault) + // "face" is a misnomer because the object describes a set of cell + // indices, but we go with the conventions of the parser here... + for (size_t cartElemIdx : face) + cartElemFaultIdx_[cartElemIdx] = faultIdx; + } +} + +template +std::vector +EclGenericThresholdPressure:: +getRestartVector() const +{ + if (!enableThresholdPressure_) + return {}; + + std::vector result(numEquilRegions_ * numEquilRegions_, 0.0); + const auto& simConfig = eclState_.getSimulationConfig(); + const auto& thpres = simConfig.getThresholdPressure(); + + std::size_t idx = 0; + for (unsigned j = 1; j <= numEquilRegions_; ++j) { + for (unsigned i = 1; i <= numEquilRegions_; ++i, ++idx) { + if (thpres.hasRegionBarrier(i, j)) { + if (thpres.hasThresholdPressure(i, j)) { + result[idx] = thpres.getThresholdPressure(i, j); + } else { + result[idx] = this->thpresDefault_[idx]; + } + } + } + } + + return result; +} + +template +void +EclGenericThresholdPressure:: +logPressures() +{ + if (!enableThresholdPressure_) + return; + + auto lineFormat = [this](unsigned i, unsigned j, double val) + { + const auto& units = eclState_.getUnits(); + return fmt::format("{:4}{:>6}{:23}{:>6}{:24}{:>11.07}{:7}{}\n", + " ", i, + " ", j, + " ", units.from_si(UnitSystem::measure::pressure, val), + " ", units.name(UnitSystem::measure::pressure)); + }; + auto lineFormatS = [](auto s1, auto s2, auto s3) + { + return fmt::format("{:4}{:^16}{:13}{:^9}{:21}{:^18}\n", + " ", s1, " ", s2, " ", s3); + }; + + std::string str = "\nLIST OF ALL NON-ZERO THRESHOLD PRESSURES\n" + "----------------------------------------\n" + "\n"; + str += lineFormatS("FLOW FROM REGION", "TO REGION", "THRESHOLD PRESSURE"); + str += lineFormatS(std::string(16, '-'), std::string(9, '-'), std::string(18, '-')); + const auto& simConfig = eclState_.getSimulationConfig(); + const auto& thpres = simConfig.getThresholdPressure(); + + for (unsigned i = 1; i <= numEquilRegions_; ++i) { + for (unsigned j = (thpres.irreversible() ? 1 : i); j <= numEquilRegions_; ++j) { + if (thpres.hasRegionBarrier(i, j)) { + if (thpres.hasThresholdPressure(i, j)) { + str += lineFormat(i, j, thpres.getThresholdPressure(j, i)); + } else { + std::size_t idx = (j - 1) * numEquilRegions_ + (i - 1); + str += lineFormat(i, j, this->thpresDefault_[idx]); + } + } + } + } + str += lineFormatS(std::string(16, '-'), std::string(9, '-'), std::string(18, '-')); + OpmLog::note(str); +} + +} // namespace Opm +#endif diff --git a/ebos/eclgenerictracermodel.cc b/ebos/eclgenerictracermodel.cc index 290d5d038..d83b442b4 100644 --- a/ebos/eclgenerictracermodel.cc +++ b/ebos/eclgenerictracermodel.cc @@ -23,7 +23,14 @@ #include #include "eclgenerictracermodel_impl.hh" +#if HAVE_DUNE_FEM +#include +#include +#include +#endif // HAVE_DUNE_FEM + namespace Opm { + #if HAVE_DUNE_FEM template class EclGenericTracerModel>>, @@ -39,54 +46,12 @@ template class EclGenericTracerModel >, false, false>, double>; -#if HAVE_DUNE_ALUGRID - -#if HAVE_MPI - using ALUGrid3CN = Dune::ALUGrid<3, 3, Dune::cube, Dune::nonconforming, Dune::ALUGridMPIComm>; -#else - using ALUGrid3CN = Dune::ALUGrid<3, 3, Dune::cube, Dune::nonconforming, Dune::ALUGridNoComm>; -#endif //HAVE_MPI - -template class EclGenericTracerModel>>, Dune::MultipleCodimMultipleGeomTypeMapper>>>, Opm::EcfvStencil>>,false,false>, - double>; - -template class EclGenericTracerModel >, - Dune::MultipleCodimMultipleGeomTypeMapper< - Dune::Fem::GridPart2GridViewImpl< -Dune::Fem::AdaptiveLeafGridPart > >, - Opm::EcfvStencil >, - false, false>, - double>; -#endif //HAVE_DUNE_ALUGRID #else template class EclGenericTracerModel>, Dune::MultipleCodimMultipleGeomTypeMapper>>, Opm::EcfvStencil>,false,false>, double>; -#if HAVE_DUNE_ALUGRID -#if HAVE_MPI - using ALUGrid3CN = Dune::ALUGrid<3, 3, Dune::cube, Dune::nonconforming, Dune::ALUGridMPIComm>; -#else - using ALUGrid3CN = Dune::ALUGrid<3, 3, Dune::cube, Dune::nonconforming, Dune::ALUGridNoComm>; -#endif //HAVE_MPI - -template class EclGenericTracerModel>, - Dune::MultipleCodimMultipleGeomTypeMapper>>, - Opm::EcfvStencil>,false,false>, - double>; -#endif //HAVE_DUNE_ALUGRID #endif //HAVE_DUNE_FEM -template class EclGenericTracerModel, - Dune::GridView>, - Dune::MultipleCodimMultipleGeomTypeMapper>>, - Opm::EcfvStencil>,false,false>, - double>; - } // namespace Opm diff --git a/ebos/eclgenerictracermodel_impl.hh b/ebos/eclgenerictracermodel_impl.hh index be71ba663..798ba1484 100644 --- a/ebos/eclgenerictracermodel_impl.hh +++ b/ebos/eclgenerictracermodel_impl.hh @@ -39,7 +39,6 @@ #include #include -#include #include #include @@ -53,18 +52,6 @@ #include #include -#if HAVE_DUNE_FEM -#include -#include -#include -#endif // HAVE_DUNE_FEM - -#if HAVE_DUNE_ALUGRID -#include -#include -#include "alucartesianindexmapper.hh" -#endif // HAVE_DUNE_ALUGRID - #include #include diff --git a/ebos/eclgenericwriter.cc b/ebos/eclgenericwriter.cc index f85f31c43..2ade8cc99 100644 --- a/ebos/eclgenericwriter.cc +++ b/ebos/eclgenericwriter.cc @@ -23,40 +23,12 @@ #include -#include - +// need to include these before impl because of specializations. #include #include -#include -#include -#if HAVE_DUNE_ALUGRID -#include "eclalugridvanguard.hh" -#include -#include -#endif // HAVE_DUNE_ALUGRID - -#include -#include -#include - -#include -#include - -#include -#include -#include -#include -#include -#include - -#include - -#include - -#if HAVE_MPI -#include -#endif +#include +#include #if HAVE_DUNE_FEM #include @@ -64,543 +36,8 @@ #include #endif // HAVE_DUNE_FEM -#if HAVE_MPI -#include -#endif - -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include - -namespace { - -/*! - * \brief Detect whether two cells are direct vertical neighbours. - * - * I.e. have the same i and j index and all cartesian cells between them - * along the vertical column are inactive. - * - * \tparam CM The type of the cartesian index mapper. - * \param cartMapper The mapper onto cartesian indices. - * \param cartesianToActive The mapping of cartesian indices to active indices. - * \param smallGlobalIndex The cartesian cell index of the cell with smaller index - * \param largeGlobalIndex The cartesian cell index of the cell with larger index - * \return True if the cells have the same i and j indices and all cartesian cells - * between them are inactive. - */ -bool directVerticalNeighbors(const std::array& cartDims, - const std::unordered_map& cartesianToActive, - int smallGlobalIndex, int largeGlobalIndex) -{ - assert(smallGlobalIndex <= largeGlobalIndex); - std::array ijk1, ijk2; - auto globalToIjk = [cartDims](int gc) { - std::array ijk; - ijk[0] = gc % cartDims[0]; - gc /= cartDims[0]; - ijk[1] = gc % cartDims[1]; - ijk[2] = gc / cartDims[1]; - return ijk; - }; - ijk1 = globalToIjk(smallGlobalIndex); - ijk2 = globalToIjk(largeGlobalIndex); - assert(ijk2[2]>=ijk1[2]); - - if ( ijk1[0] == ijk2[0] && ijk1[1] == ijk2[1] && (ijk2[2] - ijk1[2]) > 1) - { - assert((largeGlobalIndex-smallGlobalIndex)%(cartDims[0]*cartDims[1])==0); - for ( int gi = smallGlobalIndex + cartDims[0] * cartDims[1]; gi < largeGlobalIndex; - gi += cartDims[0] * cartDims[1] ) - { - if ( cartesianToActive.find( gi ) != cartesianToActive.end() ) - { - return false; - } - } - return true; - } else - return false; -} - -std::unordered_map -getInterRegFlowsAsMap(const Opm::EclInterRegFlowMap& map) -{ - auto maps = std::unordered_map{}; - - const auto& regionNames = map.names(); - auto flows = map.getInterRegFlows(); - const auto nmap = regionNames.size(); - - maps.reserve(nmap); - for (auto mapID = 0*nmap; mapID < nmap; ++mapID) { - maps.emplace(regionNames[mapID], std::move(flows[mapID])); - } - - return maps; -} - -struct EclWriteTasklet : public Opm::TaskletInterface -{ - Opm::Action::State actionState_; - Opm::WellTestState wtestState_; - Opm::SummaryState summaryState_; - Opm::UDQState udqState_; - Opm::EclipseIO& eclIO_; - int reportStepNum_; - bool isSubStep_; - double secondsElapsed_; - Opm::RestartValue restartValue_; - bool writeDoublePrecision_; - - explicit EclWriteTasklet(const Opm::Action::State& actionState, - const Opm::WellTestState& wtestState, - const Opm::SummaryState& summaryState, - const Opm::UDQState& udqState, - Opm::EclipseIO& eclIO, - int reportStepNum, - bool isSubStep, - double secondsElapsed, - Opm::RestartValue restartValue, - bool writeDoublePrecision) - : actionState_(actionState) - , wtestState_(wtestState) - , summaryState_(summaryState) - , udqState_(udqState) - , eclIO_(eclIO) - , reportStepNum_(reportStepNum) - , isSubStep_(isSubStep) - , secondsElapsed_(secondsElapsed) - , restartValue_(std::move(restartValue)) - , writeDoublePrecision_(writeDoublePrecision) - {} - - // callback to eclIO serial writeTimeStep method - void run() - { - this->eclIO_.writeTimeStep(this->actionState_, - this->wtestState_, - this->summaryState_, - this->udqState_, - this->reportStepNum_, - this->isSubStep_, - this->secondsElapsed_, - std::move(this->restartValue_), - this->writeDoublePrecision_); - } -}; - -} - namespace Opm { -template -EclGenericWriter:: -EclGenericWriter(const Schedule& schedule, - const EclipseState& eclState, - const SummaryConfig& summaryConfig, - const Grid& grid, - const EquilGrid* equilGrid, - const GridView& gridView, - const Dune::CartesianIndexMapper& cartMapper, - const Dune::CartesianIndexMapper* equilCartMapper, - bool enableAsyncOutput, - bool enableEsmry ) - : collectToIORank_(grid, - equilGrid, - gridView, - cartMapper, - equilCartMapper, - summaryConfig.fip_regions_interreg_flow()) - , grid_ (grid) - , gridView_ (gridView) - , schedule_ (schedule) - , eclState_ (eclState) - , summaryConfig_ (summaryConfig) - , cartMapper_ (cartMapper) - , equilCartMapper_(equilCartMapper) - , equilGrid_ (equilGrid) -{ - if (this->collectToIORank_.isIORank()) { - this->eclIO_ = std::make_unique - (this->eclState_, - UgGridHelpers::createEclipseGrid(*equilGrid, eclState_.getInputGrid()), - this->schedule_, this->summaryConfig_, "", enableEsmry); - } - - // create output thread if enabled and rank is I/O rank - // async output is enabled by default if pthread are enabled - int numWorkerThreads = 0; - if (enableAsyncOutput && collectToIORank_.isIORank()) { - numWorkerThreads = 1; - } - - this->taskletRunner_.reset(new TaskletRunner(numWorkerThreads)); -} - -template -const EclipseIO& EclGenericWriter:: -eclIO() const -{ - assert(eclIO_); - return *eclIO_; -} - -template -void EclGenericWriter:: -writeInit(const std::function& map) -{ - if (collectToIORank_.isIORank()) { - std::map> integerVectors; - if (collectToIORank_.isParallel()) { - integerVectors.emplace("MPI_RANK", collectToIORank_.globalRanks()); - } - - auto cartMap = cartesianToCompressed(equilGrid_->size(0), UgGridHelpers::globalCell(*equilGrid_)); - - eclIO_->writeInitial(computeTrans_(cartMap, map), - integerVectors, - exportNncStructure_(cartMap, map)); - } - -#if HAVE_MPI - if (collectToIORank_.isParallel()) { - const auto& comm = grid_.comm(); - Opm::EclMpiSerializer ser(comm); - ser.broadcast(outputNnc_); - } -#endif -} - -template -data::Solution EclGenericWriter:: -computeTrans_(const std::unordered_map& cartesianToActive, - const std::function& map) const -{ - const auto& cartMapper = *equilCartMapper_; - const auto& cartDims = cartMapper.cartesianDimensions(); - const int globalSize = cartDims[0] * cartDims[1] * cartDims[2]; - - auto tranx = data::CellData { - UnitSystem::measure::transmissibility, - std::vector(globalSize, 0.0), - data::TargetType::INIT - }; - - auto trany = tranx; - auto tranz = tranx; - - using GlobalGridView = typename EquilGrid::LeafGridView; - const GlobalGridView& globalGridView = equilGrid_->leafGridView(); - using GlobElementMapper = Dune::MultipleCodimMultipleGeomTypeMapper; - GlobElementMapper globalElemMapper(globalGridView, Dune::mcmgElementLayout()); - - for (const auto& elem : elements(globalGridView)) { - for (const auto& is : intersections(globalGridView, elem)) { - if (!is.neighbor()) - continue; // intersection is on the domain boundary - - unsigned c1 = globalElemMapper.index(is.inside()); - unsigned c2 = globalElemMapper.index(is.outside()); - - if (c1 > c2) - continue; // we only need to handle each connection once, thank you. - - // Ordering of compressed and uncompressed index should be the same - const int cartIdx1 = cartMapper.cartesianIndex( c1 ); - const int cartIdx2 = cartMapper.cartesianIndex( c2 ); - // Ordering of compressed and uncompressed index should be the same - assert(cartIdx1 <= cartIdx2); - int gc1 = std::min(cartIdx1, cartIdx2); - int gc2 = std::max(cartIdx1, cartIdx2); - - // Re-ordering in case of non-empty mapping between equilGrid to grid - if (map) { - c1 = map(c1); // equilGridToGrid map - c2 = map(c2); - } - - if (gc2 - gc1 == 1 && cartDims[0] > 1 ) { - tranx.data[gc1] = globalTrans().transmissibility(c1, c2); - continue; // skip other if clauses as they are false, last one needs some computation - } - - if (gc2 - gc1 == cartDims[0] && cartDims[1] > 1) { - trany.data[gc1] = globalTrans().transmissibility(c1, c2); - continue; // skipt next if clause as it needs some computation - } - - if ( gc2 - gc1 == cartDims[0]*cartDims[1] || - directVerticalNeighbors(cartDims, cartesianToActive, gc1, gc2)) - tranz.data[gc1] = globalTrans().transmissibility(c1, c2); - } - } - - return { - {"TRANX", tranx}, - {"TRANY", trany}, - {"TRANZ", tranz}, - }; -} - -template -std::vector EclGenericWriter:: -exportNncStructure_(const std::unordered_map& cartesianToActive, const std::function& map) const -{ - std::size_t nx = eclState_.getInputGrid().getNX(); - std::size_t ny = eclState_.getInputGrid().getNY(); - auto nncData = eclState_.getInputNNC().input(); - const auto& unitSystem = eclState_.getDeckUnitSystem(); - - for( const auto& entry : nncData ) { - // test whether NNC is not a neighboring connection - // cell2>=cell1 holds due to sortNncAndApplyEditnnc - assert( entry.cell2 >= entry.cell1 ); - auto cellDiff = entry.cell2 - entry.cell1; - - if (cellDiff != 1 && cellDiff != nx && cellDiff != nx*ny) { - auto tt = unitSystem.from_si(UnitSystem::measure::transmissibility, entry.trans); - // Eclipse ignores NNCs (with EDITNNC applied) that are small. Seems like the threshold is 1.0e-6 - if ( tt >= 1.0e-6 ) - outputNnc_.emplace_back(entry.cell1, entry.cell2, entry.trans); - } - } - - using GlobalGridView = typename EquilGrid::LeafGridView; - const GlobalGridView& globalGridView = equilGrid_->leafGridView(); - using GlobElementMapper = Dune::MultipleCodimMultipleGeomTypeMapper; - GlobElementMapper globalElemMapper(globalGridView, Dune::mcmgElementLayout()); - - // Cartesian index mapper for the serial I/O grid - const auto& equilCartMapper = *equilCartMapper_; - const auto& cartDims = cartMapper_.cartesianDimensions(); - for (const auto& elem : elements(globalGridView)) { - for (const auto& is : intersections(globalGridView, elem)) { - if (!is.neighbor()) - continue; // intersection is on the domain boundary - - unsigned c1 = globalElemMapper.index(is.inside()); - unsigned c2 = globalElemMapper.index(is.outside()); - - if (c1 > c2) - continue; // we only need to handle each connection once, thank you. - - std::size_t cc1 = equilCartMapper.cartesianIndex( c1 ); - std::size_t cc2 = equilCartMapper.cartesianIndex( c2 ); - - if ( cc2 < cc1 ) - std::swap(cc1, cc2); - - auto cellDiff = cc2 - cc1; - - // Re-ordering in case of non-empty mapping between equilGrid to grid - if (map) { - c1 = map(c1); // equilGridToGrid map - c2 = map(c2); - } - - if (cellDiff != 1 && - cellDiff != nx && - cellDiff != nx*ny && - !directVerticalNeighbors(cartDims, cartesianToActive, cc1, cc2)) { - // We need to check whether an NNC for this face was also specified - // via the NNC keyword in the deck (i.e. in the first origNncSize entries. - auto t = globalTrans().transmissibility(c1, c2); - auto candidate = std::lower_bound(nncData.begin(), nncData.end(), NNCdata(cc1, cc2, 0.0)); - - while ( candidate != nncData.end() && candidate->cell1 == cc1 - && candidate->cell2 == cc2) { - t -= candidate->trans; - ++candidate; - } - // eclipse ignores NNCs with zero transmissibility (different threshold than for NNC - // with corresponding EDITNNC above). In addition we do set small transmissibilties - // to zero when setting up the simulator. These will be ignored here, too. - auto tt = unitSystem.from_si(UnitSystem::measure::transmissibility, std::abs(t)); - if ( tt > 1e-12 ) - outputNnc_.push_back({cc1, cc2, t}); - } - } - } - return outputNnc_; -} - -template -void EclGenericWriter:: -doWriteOutput(const int reportStepNum, - const bool isSubStep, - data::Solution&& localCellData, - data::Wells&& localWellData, - data::GroupAndNetworkValues&& localGroupAndNetworkData, - data::Aquifers&& localAquiferData, - WellTestState&& localWTestState, - const Action::State& actionState, - const UDQState& udqState, - const SummaryState& summaryState, - const std::vector& thresholdPressure, - Scalar curTime, - Scalar nextStepSize, - bool doublePrecision, - bool isFlowsn, - std::array, std::vector>>, 3>&& flowsn, - bool isFloresn, - std::array, std::vector>>, 3>&& floresn) -{ - const auto isParallel = this->collectToIORank_.isParallel(); - const bool needsReordering = this->collectToIORank_.doesNeedReordering(); - - RestartValue restartValue { - (isParallel || needsReordering) - ? this->collectToIORank_.globalCellData() - : std::move(localCellData), - - isParallel ? this->collectToIORank_.globalWellData() - : std::move(localWellData), - - isParallel ? this->collectToIORank_.globalGroupAndNetworkData() - : std::move(localGroupAndNetworkData), - - isParallel ? this->collectToIORank_.globalAquiferData() - : std::move(localAquiferData) - }; - - if (eclState_.getSimulationConfig().useThresholdPressure()) { - restartValue.addExtra("THRESHPR", UnitSystem::measure::pressure, - thresholdPressure); - } - - // Add suggested next timestep to extra data. - if (! isSubStep) { - restartValue.addExtra("OPMEXTRA", std::vector(1, nextStepSize)); - } - - // Add nnc flows and flores. - if (isFlowsn) { - const auto flowsn_global = isParallel ? this->collectToIORank_.globalFlowsn() : std::move(flowsn); - for (const auto& flows : flowsn_global) { - if (flows.first.empty()) - continue; - if (flows.first == "FLOGASN+") { - restartValue.addExtra(flows.first, UnitSystem::measure::gas_surface_rate, flows.second.second); - } - else { - restartValue.addExtra(flows.first, UnitSystem::measure::liquid_surface_rate, flows.second.second); - } - } - } - if (isFloresn) { - const auto floresn_global = isParallel ? this->collectToIORank_.globalFloresn() : std::move(floresn); - for (const auto& flores : floresn_global) { - if (flores.first.empty()) - continue; - restartValue.addExtra(flores.first, UnitSystem::measure::rate, flores.second.second); - } - } - - // first, create a tasklet to write the data for the current time - // step to disk - auto eclWriteTasklet = std::make_shared( - actionState, - isParallel ? this->collectToIORank_.globalWellTestState() : std::move(localWTestState), - summaryState, udqState, *this->eclIO_, - reportStepNum, isSubStep, curTime, std::move(restartValue), doublePrecision); - - // then, make sure that the previous I/O request has been completed - // and the number of incomplete tasklets does not increase between - // time steps - this->taskletRunner_->barrier(); - - // finally, start a new output writing job - this->taskletRunner_->dispatch(std::move(eclWriteTasklet)); -} - -template -void EclGenericWriter:: -evalSummary(const int reportStepNum, - const Scalar curTime, - const data::Wells& localWellData, - const data::WellBlockAveragePressures& localWBPData, - const data::GroupAndNetworkValues& localGroupAndNetworkData, - const std::map& localAquiferData, - const std::map, double>& blockData, - const std::map& miscSummaryData, - const std::map>& regionData, - const Inplace& inplace, - const Inplace& initialInPlace, - const EclInterRegFlowMap& interRegFlows, - SummaryState& summaryState, - UDQState& udqState) -{ - if (collectToIORank_.isIORank()) { - const auto& summary = eclIO_->summary(); - - const auto& wellData = this->collectToIORank_.isParallel() - ? this->collectToIORank_.globalWellData() - : localWellData; - - const auto& wbpData = this->collectToIORank_.isParallel() - ? this->collectToIORank_.globalWBPData() - : localWBPData; - - const auto& groupAndNetworkData = this->collectToIORank_.isParallel() - ? this->collectToIORank_.globalGroupAndNetworkData() - : localGroupAndNetworkData; - - const auto& aquiferData = this->collectToIORank_.isParallel() - ? this->collectToIORank_.globalAquiferData() - : localAquiferData; - - summary.eval(summaryState, - reportStepNum, - curTime, - wellData, - wbpData, - groupAndNetworkData, - miscSummaryData, - initialInPlace, - inplace, - regionData, - blockData, - aquiferData, - getInterRegFlowsAsMap(interRegFlows)); - - // Off-by-one-fun: The reportStepNum argument corresponds to the - // report step these results will be written to, whereas the - // argument to UDQ function evaluation corresponds to the report - // step we are currently on. - const auto udq_step = reportStepNum - 1; - - this->schedule_.getUDQConfig(udq_step) - .eval(udq_step, - this->schedule_, - this->schedule_.wellMatcher(udq_step), - summaryState, - udqState); - } - -#if HAVE_MPI - if (collectToIORank_.isParallel()) { - EclMpiSerializer ser(grid_.comm()); - ser.append(summaryState); - } -#endif -} - -template -const typename EclGenericWriter::TransmissibilityType& -EclGenericWriter:: -globalTrans() const -{ - assert (globalTrans_); - return *globalTrans_; -} - #if HAVE_DUNE_FEM template class EclGenericWriter>>, double>; - - -#ifdef HAVE_DUNE_ALUGRID -#if HAVE_MPI - using ALUGrid3CN = Dune::ALUGrid<3, 3, Dune::cube, Dune::nonconforming, Dune::ALUGridMPIComm>; -#else - using ALUGrid3CN = Dune::ALUGrid<3, 3, Dune::cube, Dune::nonconforming, Dune::ALUGridNoComm>; -#endif //HAVE_MPI - -template class EclGenericWriter>>, Dune::MultipleCodimMultipleGeomTypeMapper>>>, - double>; - -template class EclGenericWriter>, - Dune::MultipleCodimMultipleGeomTypeMapper< - Dune::Fem::GridPart2GridViewImpl< - Dune::Fem::AdaptiveLeafGridPart< - ALUGrid3CN, - Dune::PartitionIteratorType(4), - false>>>, - double>; -#endif // HAVE_DUNE_ALUGRID - #else // !HAVE_DUNE_FEM template class EclGenericWriter>, Dune::MultipleCodimMultipleGeomTypeMapper>>, double>; -#if HAVE_DUNE_ALUGRID -#if HAVE_MPI - using ALUGrid3CN = Dune::ALUGrid<3, 3, Dune::cube, Dune::nonconforming, Dune::ALUGridMPIComm>; -#else - using ALUGrid3CN = Dune::ALUGrid<3, 3, Dune::cube, Dune::nonconforming, Dune::ALUGridNoComm>; -#endif //HAVE_MPI -template class EclGenericWriter>, Dune::MultipleCodimMultipleGeomTypeMapper>>, - double>; +#endif // HAVE_DUNE_FEM -#endif // HAVE_DUNE_ALUGRID -#endif // !HAVE_DUNE_FEM - -template class EclGenericWriter, - Dune::PolyhedralGrid<3,3,double>, - Dune::GridView>, Dune::MultipleCodimMultipleGeomTypeMapper>>, - double>; - } // namespace Opm diff --git a/ebos/eclgenericwriter_impl.hh b/ebos/eclgenericwriter_impl.hh new file mode 100644 index 000000000..9955df1f7 --- /dev/null +++ b/ebos/eclgenericwriter_impl.hh @@ -0,0 +1,590 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/* + 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 2 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 . + + Consult the COPYING file in the top-level source directory of this + module for the precise wording of the license and the list of + copyright holders. +*/ +#ifndef EWOMS_ECL_GENERIC_WRITER_IMPL_HH +#define EWOMS_ECL_GENERIC_WRITER_IMPL_HH + +#include + +#include + +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include + +#if HAVE_MPI +#include +#endif + +#if HAVE_MPI +#include +#endif + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace { + +/*! + * \brief Detect whether two cells are direct vertical neighbours. + * + * I.e. have the same i and j index and all cartesian cells between them + * along the vertical column are inactive. + * + * \tparam CM The type of the cartesian index mapper. + * \param cartMapper The mapper onto cartesian indices. + * \param cartesianToActive The mapping of cartesian indices to active indices. + * \param smallGlobalIndex The cartesian cell index of the cell with smaller index + * \param largeGlobalIndex The cartesian cell index of the cell with larger index + * \return True if the cells have the same i and j indices and all cartesian cells + * between them are inactive. + */ +bool directVerticalNeighbors(const std::array& cartDims, + const std::unordered_map& cartesianToActive, + int smallGlobalIndex, int largeGlobalIndex) +{ + assert(smallGlobalIndex <= largeGlobalIndex); + std::array ijk1, ijk2; + auto globalToIjk = [cartDims](int gc) { + std::array ijk; + ijk[0] = gc % cartDims[0]; + gc /= cartDims[0]; + ijk[1] = gc % cartDims[1]; + ijk[2] = gc / cartDims[1]; + return ijk; + }; + ijk1 = globalToIjk(smallGlobalIndex); + ijk2 = globalToIjk(largeGlobalIndex); + assert(ijk2[2]>=ijk1[2]); + + if ( ijk1[0] == ijk2[0] && ijk1[1] == ijk2[1] && (ijk2[2] - ijk1[2]) > 1) + { + assert((largeGlobalIndex-smallGlobalIndex)%(cartDims[0]*cartDims[1])==0); + for ( int gi = smallGlobalIndex + cartDims[0] * cartDims[1]; gi < largeGlobalIndex; + gi += cartDims[0] * cartDims[1] ) + { + if ( cartesianToActive.find( gi ) != cartesianToActive.end() ) + { + return false; + } + } + return true; + } else + return false; +} + +std::unordered_map +getInterRegFlowsAsMap(const Opm::EclInterRegFlowMap& map) +{ + auto maps = std::unordered_map{}; + + const auto& regionNames = map.names(); + auto flows = map.getInterRegFlows(); + const auto nmap = regionNames.size(); + + maps.reserve(nmap); + for (auto mapID = 0*nmap; mapID < nmap; ++mapID) { + maps.emplace(regionNames[mapID], std::move(flows[mapID])); + } + + return maps; +} + +struct EclWriteTasklet : public Opm::TaskletInterface +{ + Opm::Action::State actionState_; + Opm::WellTestState wtestState_; + Opm::SummaryState summaryState_; + Opm::UDQState udqState_; + Opm::EclipseIO& eclIO_; + int reportStepNum_; + bool isSubStep_; + double secondsElapsed_; + Opm::RestartValue restartValue_; + bool writeDoublePrecision_; + + explicit EclWriteTasklet(const Opm::Action::State& actionState, + const Opm::WellTestState& wtestState, + const Opm::SummaryState& summaryState, + const Opm::UDQState& udqState, + Opm::EclipseIO& eclIO, + int reportStepNum, + bool isSubStep, + double secondsElapsed, + Opm::RestartValue restartValue, + bool writeDoublePrecision) + : actionState_(actionState) + , wtestState_(wtestState) + , summaryState_(summaryState) + , udqState_(udqState) + , eclIO_(eclIO) + , reportStepNum_(reportStepNum) + , isSubStep_(isSubStep) + , secondsElapsed_(secondsElapsed) + , restartValue_(std::move(restartValue)) + , writeDoublePrecision_(writeDoublePrecision) + {} + + // callback to eclIO serial writeTimeStep method + void run() + { + this->eclIO_.writeTimeStep(this->actionState_, + this->wtestState_, + this->summaryState_, + this->udqState_, + this->reportStepNum_, + this->isSubStep_, + this->secondsElapsed_, + std::move(this->restartValue_), + this->writeDoublePrecision_); + } +}; + +} + +namespace Opm { + +template +EclGenericWriter:: +EclGenericWriter(const Schedule& schedule, + const EclipseState& eclState, + const SummaryConfig& summaryConfig, + const Grid& grid, + const EquilGrid* equilGrid, + const GridView& gridView, + const Dune::CartesianIndexMapper& cartMapper, + const Dune::CartesianIndexMapper* equilCartMapper, + bool enableAsyncOutput, + bool enableEsmry ) + : collectToIORank_(grid, + equilGrid, + gridView, + cartMapper, + equilCartMapper, + summaryConfig.fip_regions_interreg_flow()) + , grid_ (grid) + , gridView_ (gridView) + , schedule_ (schedule) + , eclState_ (eclState) + , summaryConfig_ (summaryConfig) + , cartMapper_ (cartMapper) + , equilCartMapper_(equilCartMapper) + , equilGrid_ (equilGrid) +{ + if (this->collectToIORank_.isIORank()) { + this->eclIO_ = std::make_unique + (this->eclState_, + UgGridHelpers::createEclipseGrid(*equilGrid, eclState_.getInputGrid()), + this->schedule_, this->summaryConfig_, "", enableEsmry); + } + + // create output thread if enabled and rank is I/O rank + // async output is enabled by default if pthread are enabled + int numWorkerThreads = 0; + if (enableAsyncOutput && collectToIORank_.isIORank()) { + numWorkerThreads = 1; + } + + this->taskletRunner_.reset(new TaskletRunner(numWorkerThreads)); +} + +template +const EclipseIO& EclGenericWriter:: +eclIO() const +{ + assert(eclIO_); + return *eclIO_; +} + +template +void EclGenericWriter:: +writeInit(const std::function& map) +{ + if (collectToIORank_.isIORank()) { + std::map> integerVectors; + if (collectToIORank_.isParallel()) { + integerVectors.emplace("MPI_RANK", collectToIORank_.globalRanks()); + } + + auto cartMap = cartesianToCompressed(equilGrid_->size(0), UgGridHelpers::globalCell(*equilGrid_)); + + eclIO_->writeInitial(computeTrans_(cartMap, map), + integerVectors, + exportNncStructure_(cartMap, map)); + } + +#if HAVE_MPI + if (collectToIORank_.isParallel()) { + const auto& comm = grid_.comm(); + Opm::EclMpiSerializer ser(comm); + ser.broadcast(outputNnc_); + } +#endif +} + +template +data::Solution EclGenericWriter:: +computeTrans_(const std::unordered_map& cartesianToActive, + const std::function& map) const +{ + const auto& cartMapper = *equilCartMapper_; + const auto& cartDims = cartMapper.cartesianDimensions(); + const int globalSize = cartDims[0] * cartDims[1] * cartDims[2]; + + auto tranx = data::CellData { + UnitSystem::measure::transmissibility, + std::vector(globalSize, 0.0), + data::TargetType::INIT + }; + + auto trany = tranx; + auto tranz = tranx; + + using GlobalGridView = typename EquilGrid::LeafGridView; + const GlobalGridView& globalGridView = equilGrid_->leafGridView(); + using GlobElementMapper = Dune::MultipleCodimMultipleGeomTypeMapper; + GlobElementMapper globalElemMapper(globalGridView, Dune::mcmgElementLayout()); + + for (const auto& elem : elements(globalGridView)) { + for (const auto& is : intersections(globalGridView, elem)) { + if (!is.neighbor()) + continue; // intersection is on the domain boundary + + unsigned c1 = globalElemMapper.index(is.inside()); + unsigned c2 = globalElemMapper.index(is.outside()); + + if (c1 > c2) + continue; // we only need to handle each connection once, thank you. + + // Ordering of compressed and uncompressed index should be the same + const int cartIdx1 = cartMapper.cartesianIndex( c1 ); + const int cartIdx2 = cartMapper.cartesianIndex( c2 ); + // Ordering of compressed and uncompressed index should be the same + assert(cartIdx1 <= cartIdx2); + int gc1 = std::min(cartIdx1, cartIdx2); + int gc2 = std::max(cartIdx1, cartIdx2); + + // Re-ordering in case of non-empty mapping between equilGrid to grid + if (map) { + c1 = map(c1); // equilGridToGrid map + c2 = map(c2); + } + + if (gc2 - gc1 == 1 && cartDims[0] > 1 ) { + tranx.data[gc1] = globalTrans().transmissibility(c1, c2); + continue; // skip other if clauses as they are false, last one needs some computation + } + + if (gc2 - gc1 == cartDims[0] && cartDims[1] > 1) { + trany.data[gc1] = globalTrans().transmissibility(c1, c2); + continue; // skipt next if clause as it needs some computation + } + + if ( gc2 - gc1 == cartDims[0]*cartDims[1] || + directVerticalNeighbors(cartDims, cartesianToActive, gc1, gc2)) + tranz.data[gc1] = globalTrans().transmissibility(c1, c2); + } + } + + return { + {"TRANX", tranx}, + {"TRANY", trany}, + {"TRANZ", tranz}, + }; +} + +template +std::vector EclGenericWriter:: +exportNncStructure_(const std::unordered_map& cartesianToActive, const std::function& map) const +{ + std::size_t nx = eclState_.getInputGrid().getNX(); + std::size_t ny = eclState_.getInputGrid().getNY(); + auto nncData = eclState_.getInputNNC().input(); + const auto& unitSystem = eclState_.getDeckUnitSystem(); + + for( const auto& entry : nncData ) { + // test whether NNC is not a neighboring connection + // cell2>=cell1 holds due to sortNncAndApplyEditnnc + assert( entry.cell2 >= entry.cell1 ); + auto cellDiff = entry.cell2 - entry.cell1; + + if (cellDiff != 1 && cellDiff != nx && cellDiff != nx*ny) { + auto tt = unitSystem.from_si(UnitSystem::measure::transmissibility, entry.trans); + // Eclipse ignores NNCs (with EDITNNC applied) that are small. Seems like the threshold is 1.0e-6 + if ( tt >= 1.0e-6 ) + outputNnc_.emplace_back(entry.cell1, entry.cell2, entry.trans); + } + } + + using GlobalGridView = typename EquilGrid::LeafGridView; + const GlobalGridView& globalGridView = equilGrid_->leafGridView(); + using GlobElementMapper = Dune::MultipleCodimMultipleGeomTypeMapper; + GlobElementMapper globalElemMapper(globalGridView, Dune::mcmgElementLayout()); + + // Cartesian index mapper for the serial I/O grid + const auto& equilCartMapper = *equilCartMapper_; + const auto& cartDims = cartMapper_.cartesianDimensions(); + for (const auto& elem : elements(globalGridView)) { + for (const auto& is : intersections(globalGridView, elem)) { + if (!is.neighbor()) + continue; // intersection is on the domain boundary + + unsigned c1 = globalElemMapper.index(is.inside()); + unsigned c2 = globalElemMapper.index(is.outside()); + + if (c1 > c2) + continue; // we only need to handle each connection once, thank you. + + std::size_t cc1 = equilCartMapper.cartesianIndex( c1 ); + std::size_t cc2 = equilCartMapper.cartesianIndex( c2 ); + + if ( cc2 < cc1 ) + std::swap(cc1, cc2); + + auto cellDiff = cc2 - cc1; + + // Re-ordering in case of non-empty mapping between equilGrid to grid + if (map) { + c1 = map(c1); // equilGridToGrid map + c2 = map(c2); + } + + if (cellDiff != 1 && + cellDiff != nx && + cellDiff != nx*ny && + !directVerticalNeighbors(cartDims, cartesianToActive, cc1, cc2)) { + // We need to check whether an NNC for this face was also specified + // via the NNC keyword in the deck (i.e. in the first origNncSize entries. + auto t = globalTrans().transmissibility(c1, c2); + auto candidate = std::lower_bound(nncData.begin(), nncData.end(), NNCdata(cc1, cc2, 0.0)); + + while ( candidate != nncData.end() && candidate->cell1 == cc1 + && candidate->cell2 == cc2) { + t -= candidate->trans; + ++candidate; + } + // eclipse ignores NNCs with zero transmissibility (different threshold than for NNC + // with corresponding EDITNNC above). In addition we do set small transmissibilties + // to zero when setting up the simulator. These will be ignored here, too. + auto tt = unitSystem.from_si(UnitSystem::measure::transmissibility, std::abs(t)); + if ( tt > 1e-12 ) + outputNnc_.push_back({cc1, cc2, t}); + } + } + } + return outputNnc_; +} + +template +void EclGenericWriter:: +doWriteOutput(const int reportStepNum, + const bool isSubStep, + data::Solution&& localCellData, + data::Wells&& localWellData, + data::GroupAndNetworkValues&& localGroupAndNetworkData, + data::Aquifers&& localAquiferData, + WellTestState&& localWTestState, + const Action::State& actionState, + const UDQState& udqState, + const SummaryState& summaryState, + const std::vector& thresholdPressure, + Scalar curTime, + Scalar nextStepSize, + bool doublePrecision, + bool isFlowsn, + std::array, std::vector>>, 3>&& flowsn, + bool isFloresn, + std::array, std::vector>>, 3>&& floresn) +{ + const auto isParallel = this->collectToIORank_.isParallel(); + const bool needsReordering = this->collectToIORank_.doesNeedReordering(); + + RestartValue restartValue { + (isParallel || needsReordering) + ? this->collectToIORank_.globalCellData() + : std::move(localCellData), + + isParallel ? this->collectToIORank_.globalWellData() + : std::move(localWellData), + + isParallel ? this->collectToIORank_.globalGroupAndNetworkData() + : std::move(localGroupAndNetworkData), + + isParallel ? this->collectToIORank_.globalAquiferData() + : std::move(localAquiferData) + }; + + if (eclState_.getSimulationConfig().useThresholdPressure()) { + restartValue.addExtra("THRESHPR", UnitSystem::measure::pressure, + thresholdPressure); + } + + // Add suggested next timestep to extra data. + if (! isSubStep) { + restartValue.addExtra("OPMEXTRA", std::vector(1, nextStepSize)); + } + + // Add nnc flows and flores. + if (isFlowsn) { + const auto flowsn_global = isParallel ? this->collectToIORank_.globalFlowsn() : std::move(flowsn); + for (const auto& flows : flowsn_global) { + if (flows.first.empty()) + continue; + if (flows.first == "FLOGASN+") { + restartValue.addExtra(flows.first, UnitSystem::measure::gas_surface_rate, flows.second.second); + } + else { + restartValue.addExtra(flows.first, UnitSystem::measure::liquid_surface_rate, flows.second.second); + } + } + } + if (isFloresn) { + const auto floresn_global = isParallel ? this->collectToIORank_.globalFloresn() : std::move(floresn); + for (const auto& flores : floresn_global) { + if (flores.first.empty()) + continue; + restartValue.addExtra(flores.first, UnitSystem::measure::rate, flores.second.second); + } + } + + // first, create a tasklet to write the data for the current time + // step to disk + auto eclWriteTasklet = std::make_shared( + actionState, + isParallel ? this->collectToIORank_.globalWellTestState() : std::move(localWTestState), + summaryState, udqState, *this->eclIO_, + reportStepNum, isSubStep, curTime, std::move(restartValue), doublePrecision); + + // then, make sure that the previous I/O request has been completed + // and the number of incomplete tasklets does not increase between + // time steps + this->taskletRunner_->barrier(); + + // finally, start a new output writing job + this->taskletRunner_->dispatch(std::move(eclWriteTasklet)); +} + +template +void EclGenericWriter:: +evalSummary(const int reportStepNum, + const Scalar curTime, + const data::Wells& localWellData, + const data::WellBlockAveragePressures& localWBPData, + const data::GroupAndNetworkValues& localGroupAndNetworkData, + const std::map& localAquiferData, + const std::map, double>& blockData, + const std::map& miscSummaryData, + const std::map>& regionData, + const Inplace& inplace, + const Inplace& initialInPlace, + const EclInterRegFlowMap& interRegFlows, + SummaryState& summaryState, + UDQState& udqState) +{ + if (collectToIORank_.isIORank()) { + const auto& summary = eclIO_->summary(); + + const auto& wellData = this->collectToIORank_.isParallel() + ? this->collectToIORank_.globalWellData() + : localWellData; + + const auto& wbpData = this->collectToIORank_.isParallel() + ? this->collectToIORank_.globalWBPData() + : localWBPData; + + const auto& groupAndNetworkData = this->collectToIORank_.isParallel() + ? this->collectToIORank_.globalGroupAndNetworkData() + : localGroupAndNetworkData; + + const auto& aquiferData = this->collectToIORank_.isParallel() + ? this->collectToIORank_.globalAquiferData() + : localAquiferData; + + summary.eval(summaryState, + reportStepNum, + curTime, + wellData, + wbpData, + groupAndNetworkData, + miscSummaryData, + initialInPlace, + inplace, + regionData, + blockData, + aquiferData, + getInterRegFlowsAsMap(interRegFlows)); + + // Off-by-one-fun: The reportStepNum argument corresponds to the + // report step these results will be written to, whereas the + // argument to UDQ function evaluation corresponds to the report + // step we are currently on. + const auto udq_step = reportStepNum - 1; + + this->schedule_.getUDQConfig(udq_step) + .eval(udq_step, + this->schedule_, + this->schedule_.wellMatcher(udq_step), + summaryState, + udqState); + } + +#if HAVE_MPI + if (collectToIORank_.isParallel()) { + EclMpiSerializer ser(grid_.comm()); + ser.append(summaryState); + } +#endif +} + +template +const typename EclGenericWriter::TransmissibilityType& +EclGenericWriter:: +globalTrans() const +{ + assert (globalTrans_); + return *globalTrans_; +} + + +} // namespace Opm +#endif diff --git a/ebos/eclpolyhedralgridvanguard.hh b/ebos/eclpolyhedralgridvanguard.hh index 02f471a6d..2eb856119 100644 --- a/ebos/eclpolyhedralgridvanguard.hh +++ b/ebos/eclpolyhedralgridvanguard.hh @@ -92,12 +92,14 @@ public: using Grid = GetPropType; using EquilGrid = GetPropType; using GridView = GetPropType; + using CartesianIndexMapper = Dune::CartesianIndexMapper; + using EquilCartesianIndexMapper = Dune::CartesianIndexMapper; + static constexpr int dimension = Grid::dimension; + static constexpr int dimensionworld = Grid::dimensionworld; private: using GridPointer = Grid*; using EquilGridPointer = EquilGrid*; - using CartesianIndexMapper = Dune::CartesianIndexMapper; - using CartesianIndexMapperPointer = std::unique_ptr; public: using TransmissibilityType = EclTransmissibilitycallImplementationInit(); + // add a copy in standard vector format to fullfill new interface + const int* globalcellorg = this->grid().globalCell(); + int num_cells = this->gridView().size(0); + globalcell_.resize(num_cells); + for(int i=0; i < num_cells; ++i){ + globalcell_[i] = globalcellorg[i]; + } } ~EclPolyhedralGridVanguard() { - delete grid_; } /*! @@ -155,7 +163,8 @@ public: * (For parallel simulation runs.) */ void loadBalance() - { /* do nothing: PolyhedralGrid is not parallel! */ } + { /* do nothing: PolyhedralGrid is not parallel! */ + } /*! * \brief Returns the object which maps a global element index of the simulation grid @@ -173,6 +182,18 @@ public: const CartesianIndexMapper& equilCartesianIndexMapper() const { return *cartesianIndexMapper_; } + const std::vector& globalCell() + { + return globalcell_; + } + + unsigned int gridEquilIdxToGridIdx(unsigned int elemIndex) const { + return elemIndex; + } + + unsigned int gridIdxToEquilGridIdx(unsigned int elemIndex) const { + return elemIndex; + } /*! * \brief Free the memory occupied by the global transmissibility object. @@ -204,11 +225,17 @@ public: return this->cellCentroids_(this->cartesianIndexMapper()); } + std::vector cellPartition() const + { + // not required for this type of grid yet (only from bdaBridge??) + return {}; + } protected: void createGrids_() { - grid_ = new Grid(this->eclState().getInputGrid(), this->eclState().fieldProps().porv(true)); + grid_ = std::make_unique(this->eclState().getInputGrid(), this->eclState().fieldProps().porv(true)); cartesianIndexMapper_ = std::make_unique(*grid_); + this->updateGridView_(); this->updateCartesianToCompressedMapping_(); this->updateCellDepths_(); } @@ -220,10 +247,12 @@ protected: Simulator& simulator_; - GridPointer grid_; - CartesianIndexMapperPointer cartesianIndexMapper_; + std::unique_ptr grid_; + std::unique_ptr cartesianIndexMapper_; + //CartesianIndexMapperPointer cartesianIndexMapper_; std::unordered_set defunctWellNames_; + std::vector globalcell_; }; } // namespace Opm diff --git a/ebos/eclproblem.hh b/ebos/eclproblem.hh index cd270e6a9..b6448a507 100644 --- a/ebos/eclproblem.hh +++ b/ebos/eclproblem.hh @@ -34,6 +34,7 @@ #include #include +#include #include #include #include @@ -84,15 +85,6 @@ #include -#if USE_ALUGRID -#define DISABLE_ALUGRID_SFC_ORDERING 1 -#include -#elif USE_POLYHEDRALGRID -#include -#else -#include -#endif - #include #include #include @@ -108,19 +100,9 @@ namespace Opm::Properties { namespace TTag { -#if USE_ALUGRID -struct EclBaseProblem { - using InheritsFrom = std::tuple; -}; -#elif USE_POLYHEDRALGRID -struct EclBaseProblem { - using InheritsFrom = std::tuple; -}; -#else struct EclBaseProblem { using InheritsFrom = std::tuple; }; -#endif } // The class which deals with ECL wells @@ -207,6 +189,7 @@ struct NumPressurePointsEquil { }; + // Set the problem property template struct Problem { @@ -531,6 +514,10 @@ struct EnableTemperature { static constexpr bool value = true; }; +template +struct EnableMech { + static constexpr bool value = false; +}; // disable all extensions supported by black oil model. this should not really be // necessary but it makes things a bit more explicit template @@ -891,12 +878,9 @@ public: readThermalParameters_(); // Re-ordering in case of ALUGrid - std::function gridToEquilGrid; - #if USE_ALUGRID - gridToEquilGrid = [&simulator](unsigned int i) { + std::function gridToEquilGrid = [&simulator](unsigned int i) { return simulator.vanguard().gridIdxToEquilGridIdx(i); }; - #endif // USE_ALUGRID transmissibilities_.finishInit(gridToEquilGrid); const auto& initconfig = eclState.getInitConfig(); @@ -936,12 +920,9 @@ public: eclWriter_->setTransmissibilities(&simulator.problem().eclTransmissibilities()); // Re-ordering in case of ALUGrid - std::function equilGridToGrid; - #if USE_ALUGRID - equilGridToGrid = [&simulator](unsigned int i) { + std::function equilGridToGrid = [&simulator](unsigned int i) { return simulator.vanguard().gridEquilIdxToGridIdx(i); }; - #endif // USE_ALUGRID eclWriter_->writeInit(equilGridToGrid); } @@ -1030,12 +1011,9 @@ public: eclBroadcast(cc, eclState.getTransMult() ); // Re-ordering in case of ALUGrid - std::function equilGridToGrid; - #if USE_ALUGRID - equilGridToGrid = [&simulator](unsigned int i) { + std::function equilGridToGrid = [&simulator](unsigned int i) { return simulator.vanguard().gridEquilIdxToGridIdx(i); }; - #endif // USE_ALUGRID // re-compute all quantities which may possibly be affected. transmissibilities_.update(true, equilGridToGrid); @@ -1177,12 +1155,9 @@ public: int episodeIdx = this->episodeIndex(); // Re-ordering in case of Alugrid - std::function gridToEquilGrid; - #if USE_ALUGRID - gridToEquilGrid = [&simulator](unsigned int i) { + std::function gridToEquilGrid = [&simulator](unsigned int i) { return simulator.vanguard().gridIdxToEquilGridIdx(i); }; - #endif // USE_ALUGRID std::function transUp = [this,gridToEquilGrid](bool global) { diff --git a/ebos/ecltransmissibility.cc b/ebos/ecltransmissibility.cc index 54228b9ef..72f2fc2ee 100644 --- a/ebos/ecltransmissibility.cc +++ b/ebos/ecltransmissibility.cc @@ -23,25 +23,9 @@ #include #include - -#include -#include +#include #include -#include - -#if HAVE_DUNE_ALUGRID -#include -#include -#include "alucartesianindexmapper.hh" -#endif // HAVE_DUNE_ALUGRID - -#include -#include -#include -#include -#include -#include #if HAVE_DUNE_FEM #include @@ -49,1061 +33,8 @@ #include #endif -#include - -#include -#include -#include -#include -#include - -namespace { - -constexpr unsigned elemIdxShift = 32; // bits - -std::uint64_t isId(std::uint32_t elemIdx1, std::uint32_t elemIdx2) -{ - std::uint32_t elemAIdx = std::min(elemIdx1, elemIdx2); - std::uint64_t elemBIdx = std::max(elemIdx1, elemIdx2); - - return (elemBIdx< isIdReverse(const std::uint64_t& id) -{ - // Assigning an unsigned integer to a narrower type discards the most significant bits. - // See "The C programming language", section A.6.2. - // NOTE that the ordering of element A and B may have changed - std::uint32_t elemAIdx = id; - std::uint32_t elemBIdx = (id - elemAIdx) >> elemIdxShift; - - return std::make_pair(elemAIdx, elemBIdx); -} - -std::uint64_t directionalIsId(std::uint32_t elemIdx1, std::uint32_t elemIdx2) -{ - return (std::uint64_t(elemIdx1)< -EclTransmissibility:: -EclTransmissibility(const EclipseState& eclState, - const GridView& gridView, - const CartesianIndexMapper& cartMapper, - const Grid& grid, - std::function(int)> centroids, - bool enableEnergy, - bool enableDiffusivity) - : eclState_(eclState) - , gridView_(gridView) - , cartMapper_(cartMapper) - , grid_(grid) - , centroids_(centroids) - , enableEnergy_(enableEnergy) - , enableDiffusivity_(enableDiffusivity) -{ - const UnitSystem& unitSystem = eclState_.getDeckUnitSystem(); - transmissibilityThreshold_ = unitSystem.parse("Transmissibility").getSIScaling() * 1e-6; -} - -template -Scalar EclTransmissibility:: -transmissibility(unsigned elemIdx1, unsigned elemIdx2) const -{ - return trans_.at(isId(elemIdx1, elemIdx2)); -} - -template -Scalar EclTransmissibility:: -transmissibilityBoundary(unsigned elemIdx, unsigned boundaryFaceIdx) const -{ - return transBoundary_.at(std::make_pair(elemIdx, boundaryFaceIdx)); -} - -template -Scalar EclTransmissibility:: -thermalHalfTrans(unsigned insideElemIdx, unsigned outsideElemIdx) const -{ - return thermalHalfTrans_.at(directionalIsId(insideElemIdx, outsideElemIdx)); -} - -template -Scalar EclTransmissibility:: -thermalHalfTransBoundary(unsigned insideElemIdx, unsigned boundaryFaceIdx) const -{ - return thermalHalfTransBoundary_.at(std::make_pair(insideElemIdx, boundaryFaceIdx)); -} - -template -Scalar EclTransmissibility:: -diffusivity(unsigned elemIdx1, unsigned elemIdx2) const -{ - if (diffusivity_.empty()) - return 0.0; - - return diffusivity_.at(isId(elemIdx1, elemIdx2)); - -} - -template -void EclTransmissibility:: -update(bool global, const std::function& map) -{ - const auto& cartDims = cartMapper_.cartesianDimensions(); - auto& transMult = eclState_.getTransMult(); - const auto& comm = gridView_.comm(); - ElementMapper elemMapper(gridView_, Dune::mcmgElementLayout()); - - // get the ntg values, the ntg values are modified for the cells merged with minpv - const std::vector& ntg = eclState_.fieldProps().get_double("NTG"); - const bool updateDiffusivity = eclState_.getSimulationConfig().isDiffusive() && enableDiffusivity_; - unsigned numElements = elemMapper.size(); - - if (map) - extractPermeability_(map); - else - extractPermeability_(); - - // calculate the axis specific centroids of all elements - std::array, dimWorld> axisCentroids; - - for (unsigned dimIdx = 0; dimIdx < dimWorld; ++dimIdx) - axisCentroids[dimIdx].resize(numElements); - - for (const auto& elem : elements(gridView_)) { - unsigned elemIdx = elemMapper.index(elem); - - // compute the axis specific "centroids" used for the transmissibilities. for - // consistency with the flow simulator, we use the element centers as - // computed by opm-parser's Opm::EclipseGrid class for all axes. - std::array centroid = centroids_(elemIdx); - - for (unsigned axisIdx = 0; axisIdx < dimWorld; ++axisIdx) - for (unsigned dimIdx = 0; dimIdx < dimWorld; ++dimIdx) - axisCentroids[axisIdx][elemIdx][dimIdx] = centroid[dimIdx]; - } - - // reserving some space in the hashmap upfront saves quite a bit of time because - // resizes are costly for hashmaps and there would be quite a few of them if we - // would not have a rough idea of how large the final map will be (the rough idea - // is a conforming Cartesian grid). - trans_.clear(); - trans_.reserve(numElements*3*1.05); - - transBoundary_.clear(); - - // if energy is enabled, let's do the same for the "thermal half transmissibilities" - if (enableEnergy_) { - thermalHalfTrans_.clear(); - thermalHalfTrans_.reserve(numElements*6*1.05); - - thermalHalfTransBoundary_.clear(); - } - - // if diffusion is enabled, let's do the same for the "diffusivity" - if (updateDiffusivity) { - diffusivity_.clear(); - diffusivity_.reserve(numElements*3*1.05); - extractPorosity_(); - } - - // The MULTZ needs special case if the option is ALL - // Then the smallest multiplier is applied. - // Default is to apply the top and bottom multiplier - bool useSmallestMultiplier; - if (comm.rank() == 0) { - const auto& eclGrid = eclState_.getInputGrid(); - useSmallestMultiplier = eclGrid.getMultzOption() == PinchMode::ALL; - } - if (global && comm.size() > 1) { - comm.broadcast(&useSmallestMultiplier, 1, 0); - } - - // compute the transmissibilities for all intersections - for (const auto& elem : elements(gridView_)) { - unsigned elemIdx = elemMapper.index(elem); - - auto isIt = gridView_.ibegin(elem); - const auto& isEndIt = gridView_.iend(elem); - unsigned boundaryIsIdx = 0; - for (; isIt != isEndIt; ++ isIt) { - // store intersection, this might be costly - const auto& intersection = *isIt; - - // deal with grid boundaries - if (intersection.boundary()) { - // compute the transmissibilty for the boundary intersection - const auto& geometry = intersection.geometry(); - const auto& faceCenterInside = geometry.center(); - - auto faceAreaNormal = intersection.centerUnitOuterNormal(); - faceAreaNormal *= geometry.volume(); - - Scalar transBoundaryIs; - computeHalfTrans_(transBoundaryIs, - faceAreaNormal, - intersection.indexInInside(), - distanceVector_(faceCenterInside, - intersection.indexInInside(), - elemIdx, - axisCentroids), - permeability_[elemIdx]); - - // normally there would be two half-transmissibilities that would be - // averaged. on the grid boundary there only is the half - // transmissibility of the interior element. - transBoundary_[std::make_pair(elemIdx, boundaryIsIdx)] = transBoundaryIs; - - // for boundary intersections we also need to compute the thermal - // half transmissibilities - if (enableEnergy_) { - Scalar transBoundaryEnergyIs; - computeHalfDiffusivity_(transBoundaryEnergyIs, - faceAreaNormal, - distanceVector_(faceCenterInside, - intersection.indexInInside(), - elemIdx, - axisCentroids), - 1.0); - thermalHalfTransBoundary_[std::make_pair(elemIdx, boundaryIsIdx)] = - transBoundaryEnergyIs; - } - - ++ boundaryIsIdx; - continue; - } - - if (!intersection.neighbor()) { - // elements can be on process boundaries, i.e. they are not on the - // domain boundary yet they don't have neighbors. - ++ boundaryIsIdx; - continue; - } - - const auto& outsideElem = intersection.outside(); - unsigned outsideElemIdx = elemMapper.index(outsideElem); - - unsigned insideCartElemIdx = cartMapper_.cartesianIndex(elemIdx); - unsigned outsideCartElemIdx = cartMapper_.cartesianIndex(outsideElemIdx); - - // we only need to calculate a face's transmissibility - // once... - if (insideCartElemIdx > outsideCartElemIdx) - continue; - - // local indices of the faces of the inside and - // outside elements which contain the intersection - int insideFaceIdx = intersection.indexInInside(); - int outsideFaceIdx = intersection.indexInOutside(); - - if (insideFaceIdx == -1) { - // NNC. Set zero transmissibility, as it will be - // *added to* by applyNncToGridTrans_() later. - assert(outsideFaceIdx == -1); - trans_[isId(elemIdx, outsideElemIdx)] = 0.0; - continue; - } - - DimVector faceCenterInside; - DimVector faceCenterOutside; - DimVector faceAreaNormal; - - typename std::is_same::type isCpGrid; - computeFaceProperties(intersection, - elemIdx, - insideFaceIdx, - outsideElemIdx, - outsideFaceIdx, - faceCenterInside, - faceCenterOutside, - faceAreaNormal, - isCpGrid); - - Scalar halfTrans1; - Scalar halfTrans2; - - computeHalfTrans_(halfTrans1, - faceAreaNormal, - insideFaceIdx, - distanceVector_(faceCenterInside, - intersection.indexInInside(), - elemIdx, - axisCentroids), - permeability_[elemIdx]); - computeHalfTrans_(halfTrans2, - faceAreaNormal, - outsideFaceIdx, - distanceVector_(faceCenterOutside, - intersection.indexInOutside(), - outsideElemIdx, - axisCentroids), - permeability_[outsideElemIdx]); - - applyNtg_(halfTrans1, insideFaceIdx, elemIdx, ntg); - applyNtg_(halfTrans2, outsideFaceIdx, outsideElemIdx, ntg); - - // convert half transmissibilities to full face - // transmissibilities using the harmonic mean - Scalar trans; - if (std::abs(halfTrans1) < 1e-30 || std::abs(halfTrans2) < 1e-30) - // avoid division by zero - trans = 0.0; - else - trans = 1.0 / (1.0/halfTrans1 + 1.0/halfTrans2); - - // apply the full face transmissibility multipliers - // for the inside ... - - if (useSmallestMultiplier) - { - // Currently PINCH(4) is never queries and hence PINCH(4) == TOPBOT is assumed - // and in this branch PINCH(5) == ALL holds - applyAllZMultipliers_(trans, insideFaceIdx, outsideFaceIdx, insideCartElemIdx, - outsideCartElemIdx, transMult, cartDims, - /* pinchTop= */ false); - } - else - { - applyMultipliers_(trans, insideFaceIdx, insideCartElemIdx, transMult); - // ... and outside elements - applyMultipliers_(trans, outsideFaceIdx, outsideCartElemIdx, transMult); - } - - // apply the region multipliers (cf. the MULTREGT keyword) - FaceDir::DirEnum faceDir; - switch (insideFaceIdx) { - case 0: - case 1: - faceDir = FaceDir::XPlus; - break; - - case 2: - case 3: - faceDir = FaceDir::YPlus; - break; - - case 4: - case 5: - faceDir = FaceDir::ZPlus; - break; - - default: - throw std::logic_error("Could not determine a face direction"); - } - - trans *= transMult.getRegionMultiplier(insideCartElemIdx, - outsideCartElemIdx, - faceDir); - - trans_[isId(elemIdx, outsideElemIdx)] = trans; - - // update the "thermal half transmissibility" for the intersection - if (enableEnergy_) { - - Scalar halfDiffusivity1; - Scalar halfDiffusivity2; - - computeHalfDiffusivity_(halfDiffusivity1, - faceAreaNormal, - distanceVector_(faceCenterInside, - intersection.indexInInside(), - elemIdx, - axisCentroids), - 1.0); - computeHalfDiffusivity_(halfDiffusivity2, - faceAreaNormal, - distanceVector_(faceCenterOutside, - intersection.indexInOutside(), - outsideElemIdx, - axisCentroids), - 1.0); - //TODO Add support for multipliers - thermalHalfTrans_[directionalIsId(elemIdx, outsideElemIdx)] = halfDiffusivity1; - thermalHalfTrans_[directionalIsId(outsideElemIdx, elemIdx)] = halfDiffusivity2; - } - - // update the "diffusive half transmissibility" for the intersection - if (updateDiffusivity) { - - Scalar halfDiffusivity1; - Scalar halfDiffusivity2; - - computeHalfDiffusivity_(halfDiffusivity1, - faceAreaNormal, - distanceVector_(faceCenterInside, - intersection.indexInInside(), - elemIdx, - axisCentroids), - porosity_[elemIdx]); - computeHalfDiffusivity_(halfDiffusivity2, - faceAreaNormal, - distanceVector_(faceCenterOutside, - intersection.indexInOutside(), - outsideElemIdx, - axisCentroids), - porosity_[outsideElemIdx]); - - applyNtg_(halfDiffusivity1, insideFaceIdx, elemIdx, ntg); - applyNtg_(halfDiffusivity2, outsideFaceIdx, outsideElemIdx, ntg); - - //TODO Add support for multipliers - Scalar diffusivity; - if (std::abs(halfDiffusivity1) < 1e-30 || std::abs(halfDiffusivity2) < 1e-30) - // avoid division by zero - diffusivity = 0.0; - else - diffusivity = 1.0 / (1.0/halfDiffusivity1 + 1.0/halfDiffusivity2); - - - diffusivity_[isId(elemIdx, outsideElemIdx)] = diffusivity; - } - } - } - - // potentially overwrite and/or modify transmissibilities based on input from deck - updateFromEclState_(global); - - // Create mapping from global to local index - std::unordered_map globalToLocal; - - // loop over all elements (global grid) and store Cartesian index - for (const auto& elem : elements(grid_.leafGridView())) { - int elemIdx = elemMapper.index(elem); - int cartElemIdx = cartMapper_.cartesianIndex(elemIdx); - globalToLocal[cartElemIdx] = elemIdx; - } - applyEditNncToGridTrans_(globalToLocal); - applyNncToGridTrans_(globalToLocal); - applyEditNncrToGridTrans_(globalToLocal); - - //remove very small non-neighbouring transmissibilities - removeSmallNonCartesianTransmissibilities_(); -} - -template -void EclTransmissibility:: -extractPermeability_() -{ - unsigned numElem = gridView_.size(/*codim=*/0); - permeability_.resize(numElem); - - // read the intrinsic permeabilities from the eclState. Note that all arrays - // provided by eclState are one-per-cell of "uncompressed" grid, whereas the - // simulation grid might remove a few elements. (e.g. because it is distributed - // over several processes.) - const auto& fp = eclState_.fieldProps(); - if (fp.has_double("PERMX")) { - const std::vector& permxData = fp.get_double("PERMX"); - - std::vector permyData; - if (fp.has_double("PERMY")) - permyData = fp.get_double("PERMY"); - else - permyData = permxData; - - std::vector permzData; - if (fp.has_double("PERMZ")) - permzData = fp.get_double("PERMZ"); - else - permzData = permxData; - - for (size_t dofIdx = 0; dofIdx < numElem; ++ dofIdx) { - permeability_[dofIdx] = 0.0; - permeability_[dofIdx][0][0] = permxData[dofIdx]; - permeability_[dofIdx][1][1] = permyData[dofIdx]; - permeability_[dofIdx][2][2] = permzData[dofIdx]; - } - - // for now we don't care about non-diagonal entries - - } - else - throw std::logic_error("Can't read the intrinsic permeability from the ecl state. " - "(The PERM{X,Y,Z} keywords are missing)"); -} - -template -void EclTransmissibility:: -extractPermeability_(const std::function& map) -{ - unsigned numElem = gridView_.size(/*codim=*/0); - permeability_.resize(numElem); - - // read the intrinsic permeabilities from the eclState. Note that all arrays - // provided by eclState are one-per-cell of "uncompressed" grid, whereas the - // simulation grid might remove a few elements. (e.g. because it is distributed - // over several processes.) - const auto& fp = eclState_.fieldProps(); - if (fp.has_double("PERMX")) { - const std::vector& permxData = fp.get_double("PERMX"); - - std::vector permyData; - if (fp.has_double("PERMY")) - permyData = fp.get_double("PERMY"); - else - permyData = permxData; - - std::vector permzData; - if (fp.has_double("PERMZ")) - permzData = fp.get_double("PERMZ"); - else - permzData = permxData; - - for (size_t dofIdx = 0; dofIdx < numElem; ++ dofIdx) { - permeability_[dofIdx] = 0.0; - size_t inputDofIdx = map(dofIdx); - permeability_[dofIdx][0][0] = permxData[inputDofIdx]; - permeability_[dofIdx][1][1] = permyData[inputDofIdx]; - permeability_[dofIdx][2][2] = permzData[inputDofIdx]; - } - - // for now we don't care about non-diagonal entries - - } - else - throw std::logic_error("Can't read the intrinsic permeability from the ecl state. " - "(The PERM{X,Y,Z} keywords are missing)"); -} - -template -void EclTransmissibility:: -extractPorosity_() -{ - // read the intrinsic porosity from the eclState. Note that all arrays - // provided by eclState are one-per-cell of "uncompressed" grid, whereas the - // simulation grid might remove a few elements. (e.g. because it is distributed - // over several processes.) - const auto& fp = eclState_.fieldProps(); - if (fp.has_double("PORO")) { - porosity_ = fp.get_double("PORO"); - } - else - throw std::logic_error("Can't read the porosityfrom the ecl state. " - "(The PORO keywords are missing)"); -} - -template -void EclTransmissibility:: -removeSmallNonCartesianTransmissibilities_() -{ - const auto& cartDims = cartMapper_.cartesianDimensions(); - for (auto&& trans: trans_) { - if (trans.second < transmissibilityThreshold_) { - const auto& id = trans.first; - const auto& elements = isIdReverse(id); - int gc1 = std::min(cartMapper_.cartesianIndex(elements.first), cartMapper_.cartesianIndex(elements.second)); - int gc2 = std::max(cartMapper_.cartesianIndex(elements.first), cartMapper_.cartesianIndex(elements.second)); - - // only adjust the NNCs - if (gc2 - gc1 == 1 || gc2 - gc1 == cartDims[0] || gc2 - gc1 == cartDims[0]*cartDims[1]) - continue; - - //remove transmissibilities less than the threshold (by default 1e-6 in the deck's unit system) - trans.second = 0.0; - } - } -} - -template -void EclTransmissibility:: -applyAllZMultipliers_(Scalar& trans, - unsigned insideFaceIdx, - unsigned outsideFaceIdx, - unsigned insideCartElemIdx, - unsigned outsideCartElemIdx, - const TransMult& transMult, - const std::array& cartDims, - bool pinchTop) -{ - if (insideFaceIdx > 3) { // top or or bottom - assert(insideFaceIdx==5); // as insideCartElemIdx < outsideCartElemIdx holds for the Z column - assert(outsideCartElemIdx > insideCartElemIdx); - auto lastCartElemIdx = outsideCartElemIdx - cartDims[0]*cartDims[1]; - // Last multiplier using (Z+)*(Z-) - Scalar mult = transMult.getMultiplier(lastCartElemIdx , FaceDir::ZPlus) * - transMult.getMultiplier(outsideCartElemIdx , FaceDir::ZMinus); - - if ( !pinchTop ) - { - // pick the smallest multiplier using (Z+)*(Z-) while looking down - // the pillar until reaching the other end of the connection - for(auto cartElemIdx = insideCartElemIdx; cartElemIdx < lastCartElemIdx;) - { - auto multiplier = transMult.getMultiplier(cartElemIdx, FaceDir::ZPlus); - cartElemIdx += cartDims[0]*cartDims[1]; - multiplier *= transMult.getMultiplier(cartElemIdx, FaceDir::ZMinus); - mult = std::min(mult, multiplier); - } - } - - trans *= mult; - } - else - { - applyMultipliers_(trans, insideFaceIdx, insideCartElemIdx, transMult); - applyMultipliers_(trans, outsideFaceIdx, outsideCartElemIdx, transMult); - } -} - -template -void EclTransmissibility:: -updateFromEclState_(bool global) -{ - const FieldPropsManager* fp = - (global) ? &(eclState_.fieldProps()) : - &(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); -} - -template -std::array,3> -EclTransmissibility:: -createTransmissibilityArrays_(const std::array& is_tran) -{ - const auto& cartDims = cartMapper_.cartesianDimensions(); - ElementMapper elemMapper(gridView_, Dune::mcmgElementLayout()); - - auto numElem = 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 - for (const auto& elem : elements(gridView_)) { - for (const auto& intersection : intersections(gridView_, elem)) { - // store intersection, this might be costly - 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 EclTransmissibility:: -resetTransmissibilityFromArrays_(const std::array& is_tran, - const std::array,3>& trans) -{ - const auto& cartDims = cartMapper_.cartesianDimensions(); - ElementMapper elemMapper(gridView_, Dune::mcmgElementLayout()); - - // compute the transmissibilities for all intersections - for (const auto& elem : elements(gridView_)) { - for (const auto& intersection : intersections(gridView_, elem)) { - 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_[isID] = trans[0][c1]; - } - else if (gc2 - gc1 == cartDims[0] && cartDims[1] > 1) { - if (is_tran[1]) - // set simulator internal transmissibilities to values from inputTrany - trans_[isID] = trans[1][c1]; - } - else if (gc2 - gc1 == cartDims[0]*cartDims[1]) { - if (is_tran[2]) - // set simulator internal transmissibilities to values from inputTranz - trans_[isID] = trans[2][c1]; - } - //else.. We don't support modification of NNC at the moment. - } - } -} - -template -template -void EclTransmissibility:: -computeFaceProperties(const Intersection& intersection, - const int, - const int, - const int, - const int, - DimVector& faceCenterInside, - DimVector& faceCenterOutside, - DimVector& faceAreaNormal, - /*isCpGrid=*/std::false_type) const -{ - // default implementation for DUNE grids - const auto& geometry = intersection.geometry(); - faceCenterInside = geometry.center(); - faceCenterOutside = faceCenterInside; - - faceAreaNormal = intersection.centerUnitOuterNormal(); - faceAreaNormal *= geometry.volume(); -} - - -template -template -void EclTransmissibility:: -computeFaceProperties(const Intersection& intersection, - const int insideElemIdx, - const int insideFaceIdx, - const int outsideElemIdx, - const int outsideFaceIdx, - DimVector& faceCenterInside, - DimVector& faceCenterOutside, - DimVector& faceAreaNormal, - /*isCpGrid=*/std::true_type) const -{ - int faceIdx = intersection.id(); - faceCenterInside = grid_.faceCenterEcl(insideElemIdx, insideFaceIdx); - faceCenterOutside = grid_.faceCenterEcl(outsideElemIdx, outsideFaceIdx); - faceAreaNormal = grid_.faceAreaNormalEcl(faceIdx); -} - -template -std::tuple, std::vector> -EclTransmissibility:: -applyNncToGridTrans_(const std::unordered_map& cartesianToCompressed) -{ - // First scale NNCs with EDITNNC. - std::vector unprocessedNnc; - std::vector processedNnc; - const auto& nnc_input = eclState_.getInputNNC().input(); - if (nnc_input.empty()) - return std::make_tuple(processedNnc, unprocessedNnc); - - for (const auto& nncEntry : nnc_input) { - auto c1 = nncEntry.cell1; - auto c2 = nncEntry.cell2; - auto lowIt = cartesianToCompressed.find(c1); - auto highIt = cartesianToCompressed.find(c2); - int low = (lowIt == cartesianToCompressed.end())? -1 : lowIt->second; - int high = (highIt == cartesianToCompressed.end())? -1 : highIt->second; - - if (low > high) - std::swap(low, high); - - if (low == -1 && high == -1) - // Silently discard as it is not between active cells - continue; - - if (low == -1 || high == -1) { - // Discard the NNC if it is between active cell and inactive cell - std::ostringstream sstr; - sstr << "NNC between active and inactive cells (" - << low << " -> " << high << ") with globalcell is (" << c1 << "->" << c2 <<")"; - OpmLog::warning(sstr.str()); - continue; - } - - auto candidate = trans_.find(isId(low, high)); - - if (candidate == trans_.end()) - // This NNC is not resembled by the grid. Save it for later - // processing with local cell values - unprocessedNnc.push_back(nncEntry); - else { - // NNC is represented by the grid and might be a neighboring connection - // In this case the transmissibilty is added to the value already - // set or computed. - candidate->second += nncEntry.trans; - processedNnc.push_back(nncEntry); - } - } - return std::make_tuple(processedNnc, unprocessedNnc); -} - -template -void EclTransmissibility:: -applyEditNncToGridTrans_(const std::unordered_map& globalToLocal) -{ - const auto& input = eclState_.getInputNNC(); - applyEditNncToGridTransHelper_(globalToLocal, "EDITNNC", - input.edit(), - [&input](const NNCdata& nnc){ - return input.edit_location(nnc);}, - // Multiply transmissibility with EDITNNC value - [](double& trans, const double& rhs){ trans *= rhs;}); -} - -template -void EclTransmissibility:: -applyEditNncrToGridTrans_(const std::unordered_map& globalToLocal) -{ - const auto& input = eclState_.getInputNNC(); - applyEditNncToGridTransHelper_(globalToLocal, "EDITNNCR", - input.editr(), - [&input](const NNCdata& nnc){ - return input.editr_location(nnc);}, - // Replace Transmissibility with EDITNNCR value - [](double& trans, const double& rhs){ trans = rhs;}); -} - -template -void EclTransmissibility:: -applyEditNncToGridTransHelper_(const std::unordered_map& globalToLocal, - const std::string& keyword, - const std::vector& nncs, - const std::function& getLocation, - const std::function& apply) -{ - if (nncs.empty()) - return; - const auto& cartDims = cartMapper_.cartesianDimensions(); - - auto format_ijk = [&cartDims](std::size_t cell) -> std::string { - auto i = cell % cartDims[0]; cell /= cartDims[0]; - auto j = cell % cartDims[1]; - auto k = cell / cartDims[1]; - - return fmt::format("({},{},{})", i + 1,j + 1,k + 1); - }; - - auto print_warning = [&format_ijk, &getLocation, &keyword] (const NNCdata& nnc) { - const auto& location = getLocation( nnc ); - auto warning = fmt::format("Problem with {} keyword\n" - "In {} line {} \n" - "No NNC defined for connection {} -> {}", keyword, location.filename, - location.lineno, format_ijk(nnc.cell1), format_ijk(nnc.cell2)); - OpmLog::warning(keyword, warning); - }; - - // editNnc is supposed to only reference non-neighboring connections and not - // neighboring connections. Use all entries for scaling if there is an NNC. - // variable nnc incremented in loop body. - auto nnc = nncs.begin(); - auto end = nncs.end(); - std::size_t warning_count = 0; - while (nnc != end) { - auto c1 = nnc->cell1; - auto c2 = nnc->cell2; - auto lowIt = globalToLocal.find(c1); - auto highIt = globalToLocal.find(c2); - - if (lowIt == globalToLocal.end() || highIt == globalToLocal.end()) { - print_warning(*nnc); - ++nnc; - warning_count++; - continue; - } - - auto low = lowIt->second, high = highIt->second; - - if (low > high) - std::swap(low, high); - - auto candidate = trans_.find(isId(low, high)); - if (candidate == trans_.end()) { - print_warning(*nnc); - ++nnc; - warning_count++; - } - else { - // NNC exists - while (nnc!= end && c1==nnc->cell1 && c2==nnc->cell2) { - apply(candidate->second, nnc->trans); - ++nnc; - } - } - } - - if (warning_count > 0) { - auto warning = fmt::format("Problems with {} keyword\n" - "A total of {} connections not defined in grid", keyword, warning_count); - OpmLog::warning(warning); - } -} - -template -void EclTransmissibility:: -computeHalfTrans_(Scalar& halfTrans, - const DimVector& areaNormal, - int faceIdx, // in the reference element that contains the intersection - const DimVector& distance, - const DimMatrix& perm) const -{ - assert(faceIdx >= 0); - unsigned dimIdx = faceIdx/2; - assert(dimIdx < dimWorld); - halfTrans = perm[dimIdx][dimIdx]; - - Scalar val = 0; - for (unsigned i = 0; i < areaNormal.size(); ++i) - val += areaNormal[i]*distance[i]; - - halfTrans *= std::abs(val); - halfTrans /= distance.two_norm2(); -} - -template -void EclTransmissibility:: -computeHalfDiffusivity_(Scalar& halfDiff, - const DimVector& areaNormal, - const DimVector& distance, - const Scalar& poro) const -{ - halfDiff = poro; - Scalar val = 0; - for (unsigned i = 0; i < areaNormal.size(); ++i) - val += areaNormal[i]*distance[i]; - - halfDiff *= std::abs(val); - halfDiff /= distance.two_norm2(); -} - -template -typename EclTransmissibility::DimVector -EclTransmissibility:: -distanceVector_(const DimVector& center, - int faceIdx, // in the reference element that contains the intersection - unsigned elemIdx, - const std::array, dimWorld>& axisCentroids) const -{ - assert(faceIdx >= 0); - unsigned dimIdx = faceIdx/2; - assert(dimIdx < dimWorld); - DimVector x = center; - x -= axisCentroids[dimIdx][elemIdx]; - - return x; -} - -template -void EclTransmissibility:: -applyMultipliers_(Scalar& trans, - unsigned faceIdx, - unsigned cartElemIdx, - const TransMult& transMult) const -{ - // apply multiplyer for the transmissibility of the face. (the - // face index is the index of the reference-element face which - // contains the intersection of interest.) - switch (faceIdx) { - case 0: // left - trans *= transMult.getMultiplier(cartElemIdx, FaceDir::XMinus); - break; - case 1: // right - trans *= transMult.getMultiplier(cartElemIdx, FaceDir::XPlus); - break; - - case 2: // front - trans *= transMult.getMultiplier(cartElemIdx, FaceDir::YMinus); - break; - case 3: // back - trans *= transMult.getMultiplier(cartElemIdx, FaceDir::YPlus); - break; - - case 4: // bottom - trans *= transMult.getMultiplier(cartElemIdx, FaceDir::ZMinus); - break; - case 5: // top - trans *= transMult.getMultiplier(cartElemIdx, FaceDir::ZPlus); - break; - } -} - -template -void EclTransmissibility:: -applyNtg_(Scalar& trans, - unsigned faceIdx, - unsigned elemIdx, - const std::vector& ntg) const -{ - // apply multiplyer for the transmissibility of the face. (the - // face index is the index of the reference-element face which - // contains the intersection of interest.) - switch (faceIdx) { - case 0: // left - trans *= ntg[elemIdx]; - break; - case 1: // right - trans *= ntg[elemIdx]; - break; - - case 2: // front - trans *= ntg[elemIdx]; - break; - case 3: // back - trans *= ntg[elemIdx]; - break; - - // NTG does not apply to top and bottom faces - } -} - #ifdef HAVE_DUNE_FEM template class EclTransmissibility>>, @@ -1122,65 +53,14 @@ template class EclTransmissibility > >, - Dune::CartesianIndexMapper, + Dune::CartesianIndexMapper, double>; -#if HAVE_DUNE_ALUGRID -#if HAVE_MPI - using ALUGrid3CN = Dune::ALUGrid<3, 3, Dune::cube, Dune::nonconforming, Dune::ALUGridMPIComm>; -#else - using ALUGrid3CN = Dune::ALUGrid<3, 3, Dune::cube, Dune::nonconforming, Dune::ALUGridNoComm>; -#endif //HAVE_MPI - -template class EclTransmissibility>>, - Dune::MultipleCodimMultipleGeomTypeMapper>>>, - Dune::CartesianIndexMapper, - double>; -template class EclTransmissibility >, - Dune::MultipleCodimMultipleGeomTypeMapper< - Dune::Fem::GridPart2GridViewImpl< - Dune::Fem::AdaptiveLeafGridPart< - ALUGrid3CN, - Dune::PartitionIteratorType(4), - false> > >, - Dune::CartesianIndexMapper, - double>; -#endif //HAVE_DUNE_ALUGRID - #else // !DUNE_FEM template class EclTransmissibility>, Dune::MultipleCodimMultipleGeomTypeMapper>>, Dune::CartesianIndexMapper, double>; - -#if HAVE_DUNE_ALUGRID -#if HAVE_MPI - using ALUGrid3CN = Dune::ALUGrid<3, 3, Dune::cube, Dune::nonconforming, Dune::ALUGridMPIComm>; -#else - using ALUGrid3CN = Dune::ALUGrid<3, 3, Dune::cube, Dune::nonconforming, Dune::ALUGridNoComm>; -#endif //HAVE_MPI - -template class -EclTransmissibility - >, - Dune::MultipleCodimMultipleGeomTypeMapper< - Dune::GridView> - >, - Dune::CartesianIndexMapper, - double>; -#endif //HAVE_DUNE_ALUGRID #endif //HAVE_DUNE_FEM -template class EclTransmissibility, -Dune::GridView>, Dune::MultipleCodimMultipleGeomTypeMapper>>, -Dune::CartesianIndexMapper>, - double>; - } // namespace Opm diff --git a/ebos/ecltransmissibility_impl.hh b/ebos/ecltransmissibility_impl.hh new file mode 100644 index 000000000..2a55fb09d --- /dev/null +++ b/ebos/ecltransmissibility_impl.hh @@ -0,0 +1,1097 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/* + 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 2 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 . + + Consult the COPYING file in the top-level source directory of this + module for the precise wording of the license and the list of + copyright holders. +*/ +#ifndef EWOMS_ECL_TRANSMISSIBILITY_IMPL_HH +#define EWOMS_ECL_TRANSMISSIBILITY_IMPL_HH + +#include + +#include +#include + +#include + +#include +#include +#include +#include +#include +#include + +#include + +#include +#include +#include +#include +#include + +namespace { + +constexpr unsigned elemIdxShift = 32; // bits + +std::uint64_t isId(std::uint32_t elemIdx1, std::uint32_t elemIdx2) +{ + std::uint32_t elemAIdx = std::min(elemIdx1, elemIdx2); + std::uint64_t elemBIdx = std::max(elemIdx1, elemIdx2); + + return (elemBIdx< isIdReverse(const std::uint64_t& id) +{ + // Assigning an unsigned integer to a narrower type discards the most significant bits. + // See "The C programming language", section A.6.2. + // NOTE that the ordering of element A and B may have changed + std::uint32_t elemAIdx = id; + std::uint32_t elemBIdx = (id - elemAIdx) >> elemIdxShift; + + return std::make_pair(elemAIdx, elemBIdx); +} + +std::uint64_t directionalIsId(std::uint32_t elemIdx1, std::uint32_t elemIdx2) +{ + return (std::uint64_t(elemIdx1)< +EclTransmissibility:: +EclTransmissibility(const EclipseState& eclState, + const GridView& gridView, + const CartesianIndexMapper& cartMapper, + const Grid& grid, + std::function(int)> centroids, + bool enableEnergy, + bool enableDiffusivity) + : eclState_(eclState) + , gridView_(gridView) + , cartMapper_(cartMapper) + , grid_(grid) + , centroids_(centroids) + , enableEnergy_(enableEnergy) + , enableDiffusivity_(enableDiffusivity) +{ + const UnitSystem& unitSystem = eclState_.getDeckUnitSystem(); + transmissibilityThreshold_ = unitSystem.parse("Transmissibility").getSIScaling() * 1e-6; +} + +template +Scalar EclTransmissibility:: +transmissibility(unsigned elemIdx1, unsigned elemIdx2) const +{ + return trans_.at(isId(elemIdx1, elemIdx2)); +} + +template +Scalar EclTransmissibility:: +transmissibilityBoundary(unsigned elemIdx, unsigned boundaryFaceIdx) const +{ + return transBoundary_.at(std::make_pair(elemIdx, boundaryFaceIdx)); +} + +template +Scalar EclTransmissibility:: +thermalHalfTrans(unsigned insideElemIdx, unsigned outsideElemIdx) const +{ + return thermalHalfTrans_.at(directionalIsId(insideElemIdx, outsideElemIdx)); +} + +template +Scalar EclTransmissibility:: +thermalHalfTransBoundary(unsigned insideElemIdx, unsigned boundaryFaceIdx) const +{ + return thermalHalfTransBoundary_.at(std::make_pair(insideElemIdx, boundaryFaceIdx)); +} + +template +Scalar EclTransmissibility:: +diffusivity(unsigned elemIdx1, unsigned elemIdx2) const +{ + if (diffusivity_.empty()) + return 0.0; + + return diffusivity_.at(isId(elemIdx1, elemIdx2)); + +} + +template +void EclTransmissibility:: +update(bool global, const std::function& map) +{ + const auto& cartDims = cartMapper_.cartesianDimensions(); + auto& transMult = eclState_.getTransMult(); + const auto& comm = gridView_.comm(); + ElementMapper elemMapper(gridView_, Dune::mcmgElementLayout()); + + // get the ntg values, the ntg values are modified for the cells merged with minpv + const std::vector& ntg = eclState_.fieldProps().get_double("NTG"); + const bool updateDiffusivity = eclState_.getSimulationConfig().isDiffusive() && enableDiffusivity_; + unsigned numElements = elemMapper.size(); + + if (map) + extractPermeability_(map); + else + extractPermeability_(); + + // calculate the axis specific centroids of all elements + std::array, dimWorld> axisCentroids; + + for (unsigned dimIdx = 0; dimIdx < dimWorld; ++dimIdx) + axisCentroids[dimIdx].resize(numElements); + + for (const auto& elem : elements(gridView_)) { + unsigned elemIdx = elemMapper.index(elem); + + // compute the axis specific "centroids" used for the transmissibilities. for + // consistency with the flow simulator, we use the element centers as + // computed by opm-parser's Opm::EclipseGrid class for all axes. + std::array centroid = centroids_(elemIdx); + + for (unsigned axisIdx = 0; axisIdx < dimWorld; ++axisIdx) + for (unsigned dimIdx = 0; dimIdx < dimWorld; ++dimIdx) + axisCentroids[axisIdx][elemIdx][dimIdx] = centroid[dimIdx]; + } + + // reserving some space in the hashmap upfront saves quite a bit of time because + // resizes are costly for hashmaps and there would be quite a few of them if we + // would not have a rough idea of how large the final map will be (the rough idea + // is a conforming Cartesian grid). + trans_.clear(); + trans_.reserve(numElements*3*1.05); + + transBoundary_.clear(); + + // if energy is enabled, let's do the same for the "thermal half transmissibilities" + if (enableEnergy_) { + thermalHalfTrans_.clear(); + thermalHalfTrans_.reserve(numElements*6*1.05); + + thermalHalfTransBoundary_.clear(); + } + + // if diffusion is enabled, let's do the same for the "diffusivity" + if (updateDiffusivity) { + diffusivity_.clear(); + diffusivity_.reserve(numElements*3*1.05); + extractPorosity_(); + } + + // The MULTZ needs special case if the option is ALL + // Then the smallest multiplier is applied. + // Default is to apply the top and bottom multiplier + bool useSmallestMultiplier; + if (comm.rank() == 0) { + const auto& eclGrid = eclState_.getInputGrid(); + useSmallestMultiplier = eclGrid.getMultzOption() == PinchMode::ALL; + } + if (global && comm.size() > 1) { + comm.broadcast(&useSmallestMultiplier, 1, 0); + } + + // compute the transmissibilities for all intersections + for (const auto& elem : elements(gridView_)) { + unsigned elemIdx = elemMapper.index(elem); + + auto isIt = gridView_.ibegin(elem); + const auto& isEndIt = gridView_.iend(elem); + unsigned boundaryIsIdx = 0; + for (; isIt != isEndIt; ++ isIt) { + // store intersection, this might be costly + const auto& intersection = *isIt; + + // deal with grid boundaries + if (intersection.boundary()) { + // compute the transmissibilty for the boundary intersection + const auto& geometry = intersection.geometry(); + const auto& faceCenterInside = geometry.center(); + + auto faceAreaNormal = intersection.centerUnitOuterNormal(); + faceAreaNormal *= geometry.volume(); + + Scalar transBoundaryIs; + computeHalfTrans_(transBoundaryIs, + faceAreaNormal, + intersection.indexInInside(), + distanceVector_(faceCenterInside, + intersection.indexInInside(), + elemIdx, + axisCentroids), + permeability_[elemIdx]); + + // normally there would be two half-transmissibilities that would be + // averaged. on the grid boundary there only is the half + // transmissibility of the interior element. + transBoundary_[std::make_pair(elemIdx, boundaryIsIdx)] = transBoundaryIs; + + // for boundary intersections we also need to compute the thermal + // half transmissibilities + if (enableEnergy_) { + Scalar transBoundaryEnergyIs; + computeHalfDiffusivity_(transBoundaryEnergyIs, + faceAreaNormal, + distanceVector_(faceCenterInside, + intersection.indexInInside(), + elemIdx, + axisCentroids), + 1.0); + thermalHalfTransBoundary_[std::make_pair(elemIdx, boundaryIsIdx)] = + transBoundaryEnergyIs; + } + + ++ boundaryIsIdx; + continue; + } + + if (!intersection.neighbor()) { + // elements can be on process boundaries, i.e. they are not on the + // domain boundary yet they don't have neighbors. + ++ boundaryIsIdx; + continue; + } + + const auto& outsideElem = intersection.outside(); + unsigned outsideElemIdx = elemMapper.index(outsideElem); + + unsigned insideCartElemIdx = cartMapper_.cartesianIndex(elemIdx); + unsigned outsideCartElemIdx = cartMapper_.cartesianIndex(outsideElemIdx); + + // we only need to calculate a face's transmissibility + // once... + if (insideCartElemIdx > outsideCartElemIdx) + continue; + + // local indices of the faces of the inside and + // outside elements which contain the intersection + int insideFaceIdx = intersection.indexInInside(); + int outsideFaceIdx = intersection.indexInOutside(); + + if (insideFaceIdx == -1) { + // NNC. Set zero transmissibility, as it will be + // *added to* by applyNncToGridTrans_() later. + assert(outsideFaceIdx == -1); + trans_[isId(elemIdx, outsideElemIdx)] = 0.0; + continue; + } + + DimVector faceCenterInside; + DimVector faceCenterOutside; + DimVector faceAreaNormal; + + typename std::is_same::type isCpGrid; + computeFaceProperties(intersection, + elemIdx, + insideFaceIdx, + outsideElemIdx, + outsideFaceIdx, + faceCenterInside, + faceCenterOutside, + faceAreaNormal, + isCpGrid); + + Scalar halfTrans1; + Scalar halfTrans2; + + computeHalfTrans_(halfTrans1, + faceAreaNormal, + insideFaceIdx, + distanceVector_(faceCenterInside, + intersection.indexInInside(), + elemIdx, + axisCentroids), + permeability_[elemIdx]); + computeHalfTrans_(halfTrans2, + faceAreaNormal, + outsideFaceIdx, + distanceVector_(faceCenterOutside, + intersection.indexInOutside(), + outsideElemIdx, + axisCentroids), + permeability_[outsideElemIdx]); + + applyNtg_(halfTrans1, insideFaceIdx, elemIdx, ntg); + applyNtg_(halfTrans2, outsideFaceIdx, outsideElemIdx, ntg); + + // convert half transmissibilities to full face + // transmissibilities using the harmonic mean + Scalar trans; + if (std::abs(halfTrans1) < 1e-30 || std::abs(halfTrans2) < 1e-30) + // avoid division by zero + trans = 0.0; + else + trans = 1.0 / (1.0/halfTrans1 + 1.0/halfTrans2); + + // apply the full face transmissibility multipliers + // for the inside ... + + if (useSmallestMultiplier) + { + // Currently PINCH(4) is never queries and hence PINCH(4) == TOPBOT is assumed + // and in this branch PINCH(5) == ALL holds + applyAllZMultipliers_(trans, insideFaceIdx, outsideFaceIdx, insideCartElemIdx, + outsideCartElemIdx, transMult, cartDims, + /* pinchTop= */ false); + } + else + { + applyMultipliers_(trans, insideFaceIdx, insideCartElemIdx, transMult); + // ... and outside elements + applyMultipliers_(trans, outsideFaceIdx, outsideCartElemIdx, transMult); + } + + // apply the region multipliers (cf. the MULTREGT keyword) + FaceDir::DirEnum faceDir; + switch (insideFaceIdx) { + case 0: + case 1: + faceDir = FaceDir::XPlus; + break; + + case 2: + case 3: + faceDir = FaceDir::YPlus; + break; + + case 4: + case 5: + faceDir = FaceDir::ZPlus; + break; + + default: + throw std::logic_error("Could not determine a face direction"); + } + + trans *= transMult.getRegionMultiplier(insideCartElemIdx, + outsideCartElemIdx, + faceDir); + + trans_[isId(elemIdx, outsideElemIdx)] = trans; + + // update the "thermal half transmissibility" for the intersection + if (enableEnergy_) { + + Scalar halfDiffusivity1; + Scalar halfDiffusivity2; + + computeHalfDiffusivity_(halfDiffusivity1, + faceAreaNormal, + distanceVector_(faceCenterInside, + intersection.indexInInside(), + elemIdx, + axisCentroids), + 1.0); + computeHalfDiffusivity_(halfDiffusivity2, + faceAreaNormal, + distanceVector_(faceCenterOutside, + intersection.indexInOutside(), + outsideElemIdx, + axisCentroids), + 1.0); + //TODO Add support for multipliers + thermalHalfTrans_[directionalIsId(elemIdx, outsideElemIdx)] = halfDiffusivity1; + thermalHalfTrans_[directionalIsId(outsideElemIdx, elemIdx)] = halfDiffusivity2; + } + + // update the "diffusive half transmissibility" for the intersection + if (updateDiffusivity) { + + Scalar halfDiffusivity1; + Scalar halfDiffusivity2; + + computeHalfDiffusivity_(halfDiffusivity1, + faceAreaNormal, + distanceVector_(faceCenterInside, + intersection.indexInInside(), + elemIdx, + axisCentroids), + porosity_[elemIdx]); + computeHalfDiffusivity_(halfDiffusivity2, + faceAreaNormal, + distanceVector_(faceCenterOutside, + intersection.indexInOutside(), + outsideElemIdx, + axisCentroids), + porosity_[outsideElemIdx]); + + applyNtg_(halfDiffusivity1, insideFaceIdx, elemIdx, ntg); + applyNtg_(halfDiffusivity2, outsideFaceIdx, outsideElemIdx, ntg); + + //TODO Add support for multipliers + Scalar diffusivity; + if (std::abs(halfDiffusivity1) < 1e-30 || std::abs(halfDiffusivity2) < 1e-30) + // avoid division by zero + diffusivity = 0.0; + else + diffusivity = 1.0 / (1.0/halfDiffusivity1 + 1.0/halfDiffusivity2); + + + diffusivity_[isId(elemIdx, outsideElemIdx)] = diffusivity; + } + } + } + + // potentially overwrite and/or modify transmissibilities based on input from deck + updateFromEclState_(global); + + // Create mapping from global to local index + std::unordered_map globalToLocal; + + // loop over all elements (global grid) and store Cartesian index + for (const auto& elem : elements(grid_.leafGridView())) { + int elemIdx = elemMapper.index(elem); + int cartElemIdx = cartMapper_.cartesianIndex(elemIdx); + globalToLocal[cartElemIdx] = elemIdx; + } + applyEditNncToGridTrans_(globalToLocal); + applyNncToGridTrans_(globalToLocal); + applyEditNncrToGridTrans_(globalToLocal); + + //remove very small non-neighbouring transmissibilities + removeSmallNonCartesianTransmissibilities_(); +} + +template +void EclTransmissibility:: +extractPermeability_() +{ + unsigned numElem = gridView_.size(/*codim=*/0); + permeability_.resize(numElem); + + // read the intrinsic permeabilities from the eclState. Note that all arrays + // provided by eclState are one-per-cell of "uncompressed" grid, whereas the + // simulation grid might remove a few elements. (e.g. because it is distributed + // over several processes.) + const auto& fp = eclState_.fieldProps(); + if (fp.has_double("PERMX")) { + const std::vector& permxData = fp.get_double("PERMX"); + + std::vector permyData; + if (fp.has_double("PERMY")) + permyData = fp.get_double("PERMY"); + else + permyData = permxData; + + std::vector permzData; + if (fp.has_double("PERMZ")) + permzData = fp.get_double("PERMZ"); + else + permzData = permxData; + + for (size_t dofIdx = 0; dofIdx < numElem; ++ dofIdx) { + permeability_[dofIdx] = 0.0; + permeability_[dofIdx][0][0] = permxData[dofIdx]; + permeability_[dofIdx][1][1] = permyData[dofIdx]; + permeability_[dofIdx][2][2] = permzData[dofIdx]; + } + + // for now we don't care about non-diagonal entries + + } + else + throw std::logic_error("Can't read the intrinsic permeability from the ecl state. " + "(The PERM{X,Y,Z} keywords are missing)"); +} + +template +void EclTransmissibility:: +extractPermeability_(const std::function& map) +{ + unsigned numElem = gridView_.size(/*codim=*/0); + permeability_.resize(numElem); + + // read the intrinsic permeabilities from the eclState. Note that all arrays + // provided by eclState are one-per-cell of "uncompressed" grid, whereas the + // simulation grid might remove a few elements. (e.g. because it is distributed + // over several processes.) + const auto& fp = eclState_.fieldProps(); + if (fp.has_double("PERMX")) { + const std::vector& permxData = fp.get_double("PERMX"); + + std::vector permyData; + if (fp.has_double("PERMY")) + permyData = fp.get_double("PERMY"); + else + permyData = permxData; + + std::vector permzData; + if (fp.has_double("PERMZ")) + permzData = fp.get_double("PERMZ"); + else + permzData = permxData; + + for (size_t dofIdx = 0; dofIdx < numElem; ++ dofIdx) { + permeability_[dofIdx] = 0.0; + size_t inputDofIdx = map(dofIdx); + permeability_[dofIdx][0][0] = permxData[inputDofIdx]; + permeability_[dofIdx][1][1] = permyData[inputDofIdx]; + permeability_[dofIdx][2][2] = permzData[inputDofIdx]; + } + + // for now we don't care about non-diagonal entries + + } + else + throw std::logic_error("Can't read the intrinsic permeability from the ecl state. " + "(The PERM{X,Y,Z} keywords are missing)"); +} + +template +void EclTransmissibility:: +extractPorosity_() +{ + // read the intrinsic porosity from the eclState. Note that all arrays + // provided by eclState are one-per-cell of "uncompressed" grid, whereas the + // simulation grid might remove a few elements. (e.g. because it is distributed + // over several processes.) + const auto& fp = eclState_.fieldProps(); + if (fp.has_double("PORO")) { + porosity_ = fp.get_double("PORO"); + } + else + throw std::logic_error("Can't read the porosityfrom the ecl state. " + "(The PORO keywords are missing)"); +} + +template +void EclTransmissibility:: +removeSmallNonCartesianTransmissibilities_() +{ + const auto& cartDims = cartMapper_.cartesianDimensions(); + for (auto&& trans: trans_) { + if (trans.second < transmissibilityThreshold_) { + const auto& id = trans.first; + const auto& elements = isIdReverse(id); + int gc1 = std::min(cartMapper_.cartesianIndex(elements.first), cartMapper_.cartesianIndex(elements.second)); + int gc2 = std::max(cartMapper_.cartesianIndex(elements.first), cartMapper_.cartesianIndex(elements.second)); + + // only adjust the NNCs + if (gc2 - gc1 == 1 || gc2 - gc1 == cartDims[0] || gc2 - gc1 == cartDims[0]*cartDims[1]) + continue; + + //remove transmissibilities less than the threshold (by default 1e-6 in the deck's unit system) + trans.second = 0.0; + } + } +} + +template +void EclTransmissibility:: +applyAllZMultipliers_(Scalar& trans, + unsigned insideFaceIdx, + unsigned outsideFaceIdx, + unsigned insideCartElemIdx, + unsigned outsideCartElemIdx, + const TransMult& transMult, + const std::array& cartDims, + bool pinchTop) +{ + if (insideFaceIdx > 3) { // top or or bottom + assert(insideFaceIdx==5); // as insideCartElemIdx < outsideCartElemIdx holds for the Z column + assert(outsideCartElemIdx > insideCartElemIdx); + auto lastCartElemIdx = outsideCartElemIdx - cartDims[0]*cartDims[1]; + // Last multiplier using (Z+)*(Z-) + Scalar mult = transMult.getMultiplier(lastCartElemIdx , FaceDir::ZPlus) * + transMult.getMultiplier(outsideCartElemIdx , FaceDir::ZMinus); + + if ( !pinchTop ) + { + // pick the smallest multiplier using (Z+)*(Z-) while looking down + // the pillar until reaching the other end of the connection + for(auto cartElemIdx = insideCartElemIdx; cartElemIdx < lastCartElemIdx;) + { + auto multiplier = transMult.getMultiplier(cartElemIdx, FaceDir::ZPlus); + cartElemIdx += cartDims[0]*cartDims[1]; + multiplier *= transMult.getMultiplier(cartElemIdx, FaceDir::ZMinus); + mult = std::min(mult, multiplier); + } + } + + trans *= mult; + } + else + { + applyMultipliers_(trans, insideFaceIdx, insideCartElemIdx, transMult); + applyMultipliers_(trans, outsideFaceIdx, outsideCartElemIdx, transMult); + } +} + +template +void EclTransmissibility:: +updateFromEclState_(bool global) +{ + const FieldPropsManager* fp = + (global) ? &(eclState_.fieldProps()) : + &(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); +} + +template +std::array,3> +EclTransmissibility:: +createTransmissibilityArrays_(const std::array& is_tran) +{ + const auto& cartDims = cartMapper_.cartesianDimensions(); + ElementMapper elemMapper(gridView_, Dune::mcmgElementLayout()); + + auto numElem = 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 + for (const auto& elem : elements(gridView_)) { + for (const auto& intersection : intersections(gridView_, elem)) { + // store intersection, this might be costly + 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 EclTransmissibility:: +resetTransmissibilityFromArrays_(const std::array& is_tran, + const std::array,3>& trans) +{ + const auto& cartDims = cartMapper_.cartesianDimensions(); + ElementMapper elemMapper(gridView_, Dune::mcmgElementLayout()); + + // compute the transmissibilities for all intersections + for (const auto& elem : elements(gridView_)) { + for (const auto& intersection : intersections(gridView_, elem)) { + 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_[isID] = trans[0][c1]; + } + else if (gc2 - gc1 == cartDims[0] && cartDims[1] > 1) { + if (is_tran[1]) + // set simulator internal transmissibilities to values from inputTrany + trans_[isID] = trans[1][c1]; + } + else if (gc2 - gc1 == cartDims[0]*cartDims[1]) { + if (is_tran[2]) + // set simulator internal transmissibilities to values from inputTranz + trans_[isID] = trans[2][c1]; + } + //else.. We don't support modification of NNC at the moment. + } + } +} + +template +template +void EclTransmissibility:: +computeFaceProperties(const Intersection& intersection, + const int, + const int, + const int, + const int, + DimVector& faceCenterInside, + DimVector& faceCenterOutside, + DimVector& faceAreaNormal, + /*isCpGrid=*/std::false_type) const +{ + // default implementation for DUNE grids + const auto& geometry = intersection.geometry(); + faceCenterInside = geometry.center(); + faceCenterOutside = faceCenterInside; + + faceAreaNormal = intersection.centerUnitOuterNormal(); + faceAreaNormal *= geometry.volume(); +} + + +template +template +void EclTransmissibility:: +computeFaceProperties(const Intersection& intersection, + const int insideElemIdx, + const int insideFaceIdx, + const int outsideElemIdx, + const int outsideFaceIdx, + DimVector& faceCenterInside, + DimVector& faceCenterOutside, + DimVector& faceAreaNormal, + /*isCpGrid=*/std::true_type) const +{ + int faceIdx = intersection.id(); + faceCenterInside = grid_.faceCenterEcl(insideElemIdx, insideFaceIdx); + faceCenterOutside = grid_.faceCenterEcl(outsideElemIdx, outsideFaceIdx); + faceAreaNormal = grid_.faceAreaNormalEcl(faceIdx); +} + +template +std::tuple, std::vector> +EclTransmissibility:: +applyNncToGridTrans_(const std::unordered_map& cartesianToCompressed) +{ + // First scale NNCs with EDITNNC. + std::vector unprocessedNnc; + std::vector processedNnc; + const auto& nnc_input = eclState_.getInputNNC().input(); + if (nnc_input.empty()) + return std::make_tuple(processedNnc, unprocessedNnc); + + for (const auto& nncEntry : nnc_input) { + auto c1 = nncEntry.cell1; + auto c2 = nncEntry.cell2; + auto lowIt = cartesianToCompressed.find(c1); + auto highIt = cartesianToCompressed.find(c2); + int low = (lowIt == cartesianToCompressed.end())? -1 : lowIt->second; + int high = (highIt == cartesianToCompressed.end())? -1 : highIt->second; + + if (low > high) + std::swap(low, high); + + if (low == -1 && high == -1) + // Silently discard as it is not between active cells + continue; + + if (low == -1 || high == -1) { + // Discard the NNC if it is between active cell and inactive cell + std::ostringstream sstr; + sstr << "NNC between active and inactive cells (" + << low << " -> " << high << ") with globalcell is (" << c1 << "->" << c2 <<")"; + OpmLog::warning(sstr.str()); + continue; + } + + auto candidate = trans_.find(isId(low, high)); + + if (candidate == trans_.end()) + // This NNC is not resembled by the grid. Save it for later + // processing with local cell values + unprocessedNnc.push_back(nncEntry); + else { + // NNC is represented by the grid and might be a neighboring connection + // In this case the transmissibilty is added to the value already + // set or computed. + candidate->second += nncEntry.trans; + processedNnc.push_back(nncEntry); + } + } + return std::make_tuple(processedNnc, unprocessedNnc); +} + +template +void EclTransmissibility:: +applyEditNncToGridTrans_(const std::unordered_map& globalToLocal) +{ + const auto& input = eclState_.getInputNNC(); + applyEditNncToGridTransHelper_(globalToLocal, "EDITNNC", + input.edit(), + [&input](const NNCdata& nnc){ + return input.edit_location(nnc);}, + // Multiply transmissibility with EDITNNC value + [](double& trans, const double& rhs){ trans *= rhs;}); +} + +template +void EclTransmissibility:: +applyEditNncrToGridTrans_(const std::unordered_map& globalToLocal) +{ + const auto& input = eclState_.getInputNNC(); + applyEditNncToGridTransHelper_(globalToLocal, "EDITNNCR", + input.editr(), + [&input](const NNCdata& nnc){ + return input.editr_location(nnc);}, + // Replace Transmissibility with EDITNNCR value + [](double& trans, const double& rhs){ trans = rhs;}); +} + +template +void EclTransmissibility:: +applyEditNncToGridTransHelper_(const std::unordered_map& globalToLocal, + const std::string& keyword, + const std::vector& nncs, + const std::function& getLocation, + const std::function& apply) +{ + if (nncs.empty()) + return; + const auto& cartDims = cartMapper_.cartesianDimensions(); + + auto format_ijk = [&cartDims](std::size_t cell) -> std::string { + auto i = cell % cartDims[0]; cell /= cartDims[0]; + auto j = cell % cartDims[1]; + auto k = cell / cartDims[1]; + + return fmt::format("({},{},{})", i + 1,j + 1,k + 1); + }; + + auto print_warning = [&format_ijk, &getLocation, &keyword] (const NNCdata& nnc) { + const auto& location = getLocation( nnc ); + auto warning = fmt::format("Problem with {} keyword\n" + "In {} line {} \n" + "No NNC defined for connection {} -> {}", keyword, location.filename, + location.lineno, format_ijk(nnc.cell1), format_ijk(nnc.cell2)); + OpmLog::warning(keyword, warning); + }; + + // editNnc is supposed to only reference non-neighboring connections and not + // neighboring connections. Use all entries for scaling if there is an NNC. + // variable nnc incremented in loop body. + auto nnc = nncs.begin(); + auto end = nncs.end(); + std::size_t warning_count = 0; + while (nnc != end) { + auto c1 = nnc->cell1; + auto c2 = nnc->cell2; + auto lowIt = globalToLocal.find(c1); + auto highIt = globalToLocal.find(c2); + + if (lowIt == globalToLocal.end() || highIt == globalToLocal.end()) { + print_warning(*nnc); + ++nnc; + warning_count++; + continue; + } + + auto low = lowIt->second, high = highIt->second; + + if (low > high) + std::swap(low, high); + + auto candidate = trans_.find(isId(low, high)); + if (candidate == trans_.end()) { + print_warning(*nnc); + ++nnc; + warning_count++; + } + else { + // NNC exists + while (nnc!= end && c1==nnc->cell1 && c2==nnc->cell2) { + apply(candidate->second, nnc->trans); + ++nnc; + } + } + } + + if (warning_count > 0) { + auto warning = fmt::format("Problems with {} keyword\n" + "A total of {} connections not defined in grid", keyword, warning_count); + OpmLog::warning(warning); + } +} + +template +void EclTransmissibility:: +computeHalfTrans_(Scalar& halfTrans, + const DimVector& areaNormal, + int faceIdx, // in the reference element that contains the intersection + const DimVector& distance, + const DimMatrix& perm) const +{ + assert(faceIdx >= 0); + unsigned dimIdx = faceIdx/2; + assert(dimIdx < dimWorld); + halfTrans = perm[dimIdx][dimIdx]; + + Scalar val = 0; + for (unsigned i = 0; i < areaNormal.size(); ++i) + val += areaNormal[i]*distance[i]; + + halfTrans *= std::abs(val); + halfTrans /= distance.two_norm2(); +} + +template +void EclTransmissibility:: +computeHalfDiffusivity_(Scalar& halfDiff, + const DimVector& areaNormal, + const DimVector& distance, + const Scalar& poro) const +{ + halfDiff = poro; + Scalar val = 0; + for (unsigned i = 0; i < areaNormal.size(); ++i) + val += areaNormal[i]*distance[i]; + + halfDiff *= std::abs(val); + halfDiff /= distance.two_norm2(); +} + +template +typename EclTransmissibility::DimVector +EclTransmissibility:: +distanceVector_(const DimVector& center, + int faceIdx, // in the reference element that contains the intersection + unsigned elemIdx, + const std::array, dimWorld>& axisCentroids) const +{ + assert(faceIdx >= 0); + unsigned dimIdx = faceIdx/2; + assert(dimIdx < dimWorld); + DimVector x = center; + x -= axisCentroids[dimIdx][elemIdx]; + + return x; +} + +template +void EclTransmissibility:: +applyMultipliers_(Scalar& trans, + unsigned faceIdx, + unsigned cartElemIdx, + const TransMult& transMult) const +{ + // apply multiplyer for the transmissibility of the face. (the + // face index is the index of the reference-element face which + // contains the intersection of interest.) + switch (faceIdx) { + case 0: // left + trans *= transMult.getMultiplier(cartElemIdx, FaceDir::XMinus); + break; + case 1: // right + trans *= transMult.getMultiplier(cartElemIdx, FaceDir::XPlus); + break; + + case 2: // front + trans *= transMult.getMultiplier(cartElemIdx, FaceDir::YMinus); + break; + case 3: // back + trans *= transMult.getMultiplier(cartElemIdx, FaceDir::YPlus); + break; + + case 4: // bottom + trans *= transMult.getMultiplier(cartElemIdx, FaceDir::ZMinus); + break; + case 5: // top + trans *= transMult.getMultiplier(cartElemIdx, FaceDir::ZPlus); + break; + } +} + +template +void EclTransmissibility:: +applyNtg_(Scalar& trans, + unsigned faceIdx, + unsigned elemIdx, + const std::vector& ntg) const +{ + // apply multiplyer for the transmissibility of the face. (the + // face index is the index of the reference-element face which + // contains the intersection of interest.) + switch (faceIdx) { + case 0: // left + trans *= ntg[elemIdx]; + break; + case 1: // right + trans *= ntg[elemIdx]; + break; + + case 2: // front + trans *= ntg[elemIdx]; + break; + case 3: // back + trans *= ntg[elemIdx]; + break; + + // NTG does not apply to top and bottom faces + } +} + + +} // namespace Opm +#endif diff --git a/ebos/equil/initstateequil.cc b/ebos/equil/initstateequil.cc index dd3159b6f..53fdffa86 100644 --- a/ebos/equil/initstateequil.cc +++ b/ebos/equil/initstateequil.cc @@ -21,13 +21,22 @@ copyright holders. */ -#include "config.h" -#include "initstateequil.hh" -#include "initstateequil_impl.hh" +#include +#include +#include + +#include + +#if HAVE_DUNE_FEM +#include +#include +#include +#endif namespace Opm { namespace EQUIL { namespace DeckDependent { + #if HAVE_DUNE_FEM using GridView = Dune::Fem::GridPart2GridViewImpl< Dune::Fem::AdaptiveLeafGridPart< @@ -59,37 +68,6 @@ template InitialStateComputer, const double, const int, const bool); -#if HAVE_DUNE_ALUGRID -#if HAVE_MPI -using ALUGridComm = Dune::ALUGridMPIComm; -#else -using ALUGridComm = Dune::ALUGridNoComm; -#endif //HAVE_MPI -using ALUGrid3CN = Dune::ALUGrid<3, 3, Dune::cube, Dune::nonconforming, ALUGridComm>; -using ALUGridView = Dune::GridView>; -using ALUGridMapper = Dune::MultipleCodimMultipleGeomTypeMapper; -template class InitialStateComputer, - ALUGrid3CN, - ALUGridView, - ALUGridMapper, - Dune::CartesianIndexMapper>; - -template InitialStateComputer, - ALUGrid3CN, - ALUGridView, - ALUGridMapper, - Dune::CartesianIndexMapper>:: - InitialStateComputer(MatLaw&, - const EclipseState&, - const ALUGrid3CN&, - const ALUGridView&, - const Dune::CartesianIndexMapper&, - const double, - const int, - const bool); -#endif //HAVE_DUNE_ALUGRID - - } // namespace DeckDependent namespace Details { diff --git a/ebos/equil/initstateequil_impl.hh b/ebos/equil/initstateequil_impl.hh index 0a3a7f42f..2048f0497 100644 --- a/ebos/equil/initstateequil_impl.hh +++ b/ebos/equil/initstateequil_impl.hh @@ -29,7 +29,6 @@ #include -#include #include #include @@ -54,17 +53,6 @@ #include #include -#if HAVE_DUNE_FEM -#include -#include -#include -#endif - -#if HAVE_DUNE_ALUGRID -#include -#include -#endif // HAVE_DUNE_ALUGRID - namespace Opm { namespace EQUIL { diff --git a/flow/flow_blackoil_alugrid.cpp b/flow/flow_blackoil_alugrid.cpp index 0c2f266ed..f7a46a075 100644 --- a/flow/flow_blackoil_alugrid.cpp +++ b/flow/flow_blackoil_alugrid.cpp @@ -16,9 +16,26 @@ You should have received a copy of the GNU General Public License along with OPM. If not, see . */ -#include "config.h" +#include + +#include +#include #include +// for equilgrid in writer +// need to include this before eclgenericwriter_impl.hh due to specializations. +#include +#include + +// these are not explicitly instanced in library +#include +#include +#include +#include +#include +#include +#include + namespace Opm { namespace Properties { namespace TTag { @@ -26,6 +43,25 @@ struct EclFlowProblemAlugrid { using InheritsFrom = std::tuple; }; } + +template +struct Grid { + static const int dim = 3; +#if HAVE_MPI + using type = Dune::ALUGrid; +#else + using type = Dune::ALUGrid; +#endif +}; +// alugrid need cp grid as equilgrid +template +struct EquilGrid { + using type = Dune::CpGrid; +}; +template +struct Vanguard { + using type = Opm::EclAluGridVanguard; +}; template struct EclEnableAquifers { static constexpr bool value = false; diff --git a/flow/flow_blackoil_poly.cpp b/flow/flow_blackoil_poly.cpp deleted file mode 100644 index 8740c3af3..000000000 --- a/flow/flow_blackoil_poly.cpp +++ /dev/null @@ -1,25 +0,0 @@ -/* - Copyright 2020, NORCE AS - - 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 . -*/ -#include "config.h" -#include - -int main(int argc, char** argv) -{ - return Opm::flowEbosBlackoilPolyMain(argc, argv); -} diff --git a/flow/flow_blackoil_polyhedralgrid.cpp b/flow/flow_blackoil_polyhedralgrid.cpp new file mode 100644 index 000000000..ca1b5dd59 --- /dev/null +++ b/flow/flow_blackoil_polyhedralgrid.cpp @@ -0,0 +1,78 @@ +/* + Copyright 2020, NORCE AS + + 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 . +*/ +#include + +#include + +#include +#include +#include +#include + +// these are not explicitly instanced in library +#include +#include +#include +#include +#include +#include +#include + +namespace Opm { +namespace Properties { + namespace TTag { + struct EclFlowProblemPoly { + using InheritsFrom = std::tuple; + }; + } + + template + struct Linearizer { using type = TpfaLinearizer; }; + + template + struct LocalResidual { using type = BlackOilLocalResidualTPFA; }; + + template + struct EnableDiffusion { static constexpr bool value = false; }; + + template + struct Grid { + using type = Dune::PolyhedralGrid<3, 3>; + }; + template + struct EquilGrid { + //using type = Dune::CpGrid; + using type = GetPropType; + }; + + template + struct Vanguard { + using type = Opm::EclPolyhedralGridVanguard; + }; +} +} +int main(int argc, char** argv) +{ + using TypeTag = Opm::Properties::TTag::EclFlowProblemPoly; + auto mainObject = std::make_unique(argc, argv); + auto ret = mainObject->runStatic(); + // Destruct mainObject as the destructor calls MPI_Finalize! + mainObject.reset(); + return ret; +} diff --git a/flow/flow_ebos_blackoil_poly.cpp b/flow/flow_ebos_blackoil_poly.cpp deleted file mode 100644 index 8ddcf9a58..000000000 --- a/flow/flow_ebos_blackoil_poly.cpp +++ /dev/null @@ -1,41 +0,0 @@ -/* - Copyright 2020, NORCE AS - - 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 . -*/ -#include "config.h" -#include - -namespace Opm { -namespace Properties { - namespace TTag { - struct EclFlowProblemPoly { - using InheritsFrom = std::tuple; - }; - } -} - -int flowEbosBlackoilPolyMain(int argc, char** argv) -{ - using TypeTag = Opm::Properties::TTag::EclFlowProblemPoly; - auto mainObject = std::make_unique(argc, argv); - auto ret = mainObject->runStatic(); - // Destruct mainObject as the destructor calls MPI_Finalize! - mainObject.reset(); - return ret; -} - -} diff --git a/flow/flow_ebos_blackoil_poly.hpp b/flow/flow_ebos_blackoil_poly.hpp deleted file mode 100644 index 55f7d0f54..000000000 --- a/flow/flow_ebos_blackoil_poly.hpp +++ /dev/null @@ -1,29 +0,0 @@ -/* - Copyright 2013, 2014, 2015 SINTEF ICT, Applied Mathematics. - Copyright 2014 Dr. Blatt - HPC-Simulation-Software & Services - Copyright 2015, 2017 IRIS AS - - 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 . -*/ - -#ifndef FLOW_BLACKOIL_POLY_HPP -#define FLOW_BLACKOIL_POLY_HPP - -namespace Opm { - int flowEbosBlackoilPolyMain(int argc, char** argv); -} - -#endif diff --git a/opm/simulators/aquifers/BlackoilAquiferModel.hpp b/opm/simulators/aquifers/BlackoilAquiferModel.hpp index 37b99ae43..43c3791aa 100644 --- a/opm/simulators/aquifers/BlackoilAquiferModel.hpp +++ b/opm/simulators/aquifers/BlackoilAquiferModel.hpp @@ -37,10 +37,8 @@ #include #include -#ifdef USE_POLYHEDRALGRID #include -#endif -#if USE_ALUGRID +#if HAVE_DUNE_ALUGRID #include #endif @@ -65,14 +63,12 @@ class SupportsFaceTag {}; -#ifdef USE_POLYHEDRALGRID template<> class SupportsFaceTag> : public std::bool_constant {}; -#endif -#if USE_ALUGRID +#if HAVE_DUNE_ALUGRID template<> class SupportsFaceTag> : public std::bool_constant diff --git a/opm/simulators/linalg/ISTLSolverEbos.cpp b/opm/simulators/linalg/ISTLSolverEbos.cpp index 133c5f0c5..89f2465bf 100644 --- a/opm/simulators/linalg/ISTLSolverEbos.cpp +++ b/opm/simulators/linalg/ISTLSolverEbos.cpp @@ -43,6 +43,8 @@ #include #endif // HAVE_DUNE_ALUGRID +#include + namespace Opm { namespace detail { @@ -387,7 +389,7 @@ using CommunicationType = Dune::CollectiveCommunication; const std::vector&, \ const std::vector&, \ const size_t, const bool); - +using PolyHedralGrid3D = Dune::PolyhedralGrid<3, 3>; #if HAVE_DUNE_ALUGRID #if HAVE_MPI using ALUGrid3CN = Dune::ALUGrid<3, 3, Dune::cube, Dune::nonconforming, Dune::ALUGridMPIComm>; @@ -398,11 +400,13 @@ using CommunicationType = Dune::CollectiveCommunication; template struct BdaSolverInfo,BV>; \ INSTANCE_GRID(Dim,Dune::CpGrid) \ INSTANCE_GRID(Dim,ALUGrid3CN) \ + INSTANCE_GRID(Dim,PolyHedralGrid3D) \ INSTANCE_FLEX(Dim) #else #define INSTANCE(Dim) \ template struct BdaSolverInfo,BV>; \ INSTANCE_GRID(Dim,Dune::CpGrid) \ + INSTANCE_GRID(Dim,PolyHedralGrid3D) \ INSTANCE_FLEX(Dim) #endif #else