From e8f69223e6b2cd1169b03972f06354b49efacf22 Mon Sep 17 00:00:00 2001 From: Arne Morten Kvarving Date: Fri, 7 May 2021 08:17:22 +0200 Subject: [PATCH 1/2] changed: use constexpr value and ifs to disable code compile time instead of having runtime conditionals --- ebos/ecltransmissibility.hh | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/ebos/ecltransmissibility.hh b/ebos/ecltransmissibility.hh index 7c30ca792..4d1cd0ecf 100644 --- a/ebos/ecltransmissibility.hh +++ b/ebos/ecltransmissibility.hh @@ -71,8 +71,8 @@ class EclTransmissibility using ElementMapper = GetPropType; using Intersection = typename GridView::Intersection; - static const bool enableEnergy = getPropValue(); - static const bool enableDiffusion = getPropValue(); + static constexpr bool enableEnergy = getPropValue(); + static constexpr bool enableDiffusion = getPropValue(); // Grid and world dimension enum { dimWorld = GridView::dimensionworld }; @@ -175,7 +175,7 @@ public: transBoundary_.clear(); // if energy is enabled, let's do the same for the "thermal half transmissibilities" - if (enableEnergy) { + if constexpr (enableEnergy) { thermalHalfTrans_->clear(); thermalHalfTrans_->reserve(numElements*6*1.05); @@ -240,7 +240,7 @@ public: // for boundary intersections we also need to compute the thermal // half transmissibilities - if (enableEnergy) { + if constexpr (enableEnergy) { Scalar transBoundaryEnergyIs; computeHalfDiffusivity_(transBoundaryEnergyIs, faceAreaNormal, @@ -382,7 +382,7 @@ public: trans_[isId_(elemIdx, outsideElemIdx)] = trans; // update the "thermal half transmissibility" for the intersection - if (enableEnergy) { + if constexpr (enableEnergy) { Scalar halfDiffusivity1; Scalar halfDiffusivity2; From 7ff44d9093c4de5036e74b31d30cf3f8bfab0a88 Mon Sep 17 00:00:00 2001 From: Arne Morten Kvarving Date: Wed, 5 May 2021 12:13:25 +0200 Subject: [PATCH 2/2] ecltransmissibility: separate from typetag this allows using explicit template instantation to only compile this code per grid, not per simulator object --- CMakeLists_files.cmake | 1 + ebos/eclcpgridvanguard.hh | 18 +- ebos/eclpolyhedralgridvanguard.hh | 8 +- ebos/eclproblem.hh | 12 +- ebos/ecltransmissibility.cc | 1075 +++++++++++++++++++++++++++++ ebos/ecltransmissibility.hh | 1015 +++------------------------ ebos/eclwriter.hh | 7 +- 7 files changed, 1186 insertions(+), 950 deletions(-) create mode 100644 ebos/ecltransmissibility.cc diff --git a/CMakeLists_files.cmake b/CMakeLists_files.cmake index 75c679f8b..82983fb9e 100644 --- a/CMakeLists_files.cmake +++ b/CMakeLists_files.cmake @@ -26,6 +26,7 @@ list (APPEND MAIN_SOURCE_FILES ebos/collecttoiorank.cc ebos/eclgenericcpgridvanguard.cc ebos/eclgenericvanguard.cc + ebos/ecltransmissibility.cc opm/core/props/phaseUsageFromDeck.cpp opm/core/props/satfunc/RelpermDiagnostics.cpp opm/simulators/timestepping/SimulatorReport.cpp diff --git a/ebos/eclcpgridvanguard.hh b/ebos/eclcpgridvanguard.hh index 9f4f1e82b..4347977eb 100644 --- a/ebos/eclcpgridvanguard.hh +++ b/ebos/eclcpgridvanguard.hh @@ -27,13 +27,13 @@ #ifndef EWOMS_ECL_CP_GRID_VANGUARD_HH #define EWOMS_ECL_CP_GRID_VANGUARD_HH +#include + #include "eclbasevanguard.hh" #include "ecltransmissibility.hh" #include "femcpgridcompat.hh" #include "eclgenericcpgridvanguard.hh" -#include - namespace Opm { template class EclCpGridVanguard; @@ -89,7 +89,9 @@ public: using Grid = GetPropType; using EquilGrid = GetPropType; using GridView = GetPropType; - + using TransmissibilityType = EclTransmissibility(), + getPropValue()>; private: typedef Dune::CartesianIndexMapper CartesianIndexMapper; using Element = typename GridView::template Codim<0>::Entity; @@ -111,7 +113,7 @@ public: globalTrans_.reset(); } - const EclTransmissibility& globalTransmissibility() const + const TransmissibilityType& globalTransmissibility() const { assert( globalTrans_ != nullptr ); return *globalTrans_; @@ -152,7 +154,11 @@ protected: void allocTrans() override { - globalTrans_.reset(new EclTransmissibility(*this)); + globalTrans_.reset(new TransmissibilityType(this->eclState(), + this->gridView(), + this->cartesianIndexMapper(), + this->grid(), + this->cellCentroids())); globalTrans_->update(false); } @@ -167,7 +173,7 @@ protected: this->doFilterConnections_(this->schedule()); } - std::unique_ptr > globalTrans_; + std::unique_ptr globalTrans_; }; } // namespace Opm diff --git a/ebos/eclpolyhedralgridvanguard.hh b/ebos/eclpolyhedralgridvanguard.hh index 0657b7bcf..fcc8a1373 100644 --- a/ebos/eclpolyhedralgridvanguard.hh +++ b/ebos/eclpolyhedralgridvanguard.hh @@ -27,6 +27,8 @@ #ifndef EWOMS_ECL_POLYHEDRAL_GRID_VANGUARD_HH #define EWOMS_ECL_POLYHEDRAL_GRID_VANGUARD_HH +#include + #include "eclbasevanguard.hh" #include "ecltransmissibility.hh" @@ -76,6 +78,7 @@ class EclPolyhedralGridVanguard : public EclBaseVanguard friend class EclBaseVanguard; typedef EclBaseVanguard ParentType; + using ElementMapper = GetPropType; using Scalar = GetPropType; using Simulator = GetPropType; @@ -83,6 +86,9 @@ public: using Grid = GetPropType; using EquilGrid = GetPropType; using GridView = GetPropType; + using TransmissibilityType = EclTransmissibility(), + getPropValue()>; private: typedef Grid* GridPointer; @@ -175,7 +181,7 @@ public: std::unordered_set defunctWellNames() const { return defunctWellNames_; } - const EclTransmissibility& globalTransmissibility() const + const TransmissibilityType& globalTransmissibility() const { return simulator_.problem().eclTransmissibilities(); } diff --git a/ebos/eclproblem.hh b/ebos/eclproblem.hh index d8326c591..c168cf429 100644 --- a/ebos/eclproblem.hh +++ b/ebos/eclproblem.hh @@ -611,7 +611,7 @@ class EclProblem : public GetPropType using FluidSystem = GetPropType; using GlobalEqVector = GetPropType; using EqVector = GetPropType; - + using Vanguard = GetPropType; // Grid and world dimension enum { dim = GridView::dimension }; @@ -815,7 +815,11 @@ public: */ EclProblem(Simulator& simulator) : ParentType(simulator) - , transmissibilities_(simulator.vanguard()) + , transmissibilities_(simulator.vanguard().eclState(), + simulator.vanguard().gridView(), + simulator.vanguard().cartesianIndexMapper(), + simulator.vanguard().grid(), + simulator.vanguard().cellCentroids()) , thresholdPressures_(simulator) , wellModel_(simulator) , aquiferModel_(simulator) @@ -1514,7 +1518,7 @@ public: /*! * \brief Return a reference to the object that handles the "raw" transmissibilities. */ - const EclTransmissibility& eclTransmissibilities() const + const typename Vanguard::TransmissibilityType& eclTransmissibilities() const { return transmissibilities_; } /*! @@ -3450,7 +3454,7 @@ private: static std::string briefDescription_; std::array, 2> referencePorosity_; - EclTransmissibility transmissibilities_; + typename Vanguard::TransmissibilityType transmissibilities_; std::shared_ptr materialLawManager_; std::shared_ptr thermalLawManager_; diff --git a/ebos/ecltransmissibility.cc b/ebos/ecltransmissibility.cc new file mode 100644 index 000000000..2fec9fbf1 --- /dev/null +++ b/ebos/ecltransmissibility.cc @@ -0,0 +1,1075 @@ +// -*- 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. +*/ + +#include +#include + +#include +#include + +#include +#include + +#include +#include +#include +#include +#include + +#if HAVE_DUNE_FEM +#include +#include +#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 Dune::CartesianIndexMapper& cartMapper, + const Grid& grid, + const std::vector& centroids) + : eclState_(eclState) + , gridView_(gridView) + , cartMapper_(cartMapper) + , grid_(grid) + , centroids_(centroids) +{ + 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 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() && enableDiffusion; + unsigned numElements = elemMapper.size(); + + extractPermeability_(); + + // calculate the axis specific centroids of all elements + std::array, dimWorld> axisCentroids; + + for (unsigned dimIdx = 0; dimIdx < dimWorld; ++dimIdx) + axisCentroids[dimIdx].resize(numElements); + + auto elemIt = gridView_.template begin(); + const auto& elemEndIt = gridView_.template end(); + size_t centroidIdx = 0; + for (; elemIt != elemEndIt; ++elemIt, ++centroidIdx) { + const auto& elem = *elemIt; + 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; + if (gridView_.comm().rank() == 0) { + const auto& eclGrid = eclState_.getInputGrid(); + unsigned cartesianCellIdx = cartMapper_.cartesianIndex(elemIdx); + centroid = eclGrid.getCellCenter(cartesianCellIdx); + } else + std::copy(centroids_.begin() + centroidIdx * dimWorld, + centroids_.begin() + (centroidIdx + 1) * dimWorld, + centroid.begin()); + + 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 constexpr (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() == Opm::PinchMode::ModeEnum::ALL; + } + if (global && comm.size() > 1) { + comm.broadcast(&useSmallestMultiplier, 1, 0); + } + + // compute the transmissibilities for all intersections + elemIt = gridView_.template begin(); + for (; elemIt != elemEndIt; ++elemIt) { + const auto& elem = *elemIt; + 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 constexpr (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 constexpr (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 + const size_t cartesianSize = cartMapper_.cartesianSize(); + // reserve memory + std::vector globalToLocal(cartesianSize, -1); + + // loop over all elements (global grid) and store Cartesian index + elemIt = grid_.leafGridView().template begin<0>(); + + for (; elemIt != elemEndIt; ++elemIt) { + int elemIdx = elemMapper.index(*elemIt); + int cartElemIdx = cartMapper_.cartesianIndex(elemIdx); + globalToLocal[cartElemIdx] = elemIdx; + } + applyEditNncToGridTrans_(globalToLocal); + applyNncToGridTrans_(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:: +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 + Scalar mult = transMult.getMultiplier(lastCartElemIdx , FaceDir::ZPlus); + + if ( !pinchTop ) + { + // pick the smallest multiplier for Z+ while looking down the pillar until reaching the other end of the connection + // While Z- is not all used here. + for(auto cartElemIdx = insideCartElemIdx; cartElemIdx < lastCartElemIdx; + cartElemIdx += cartDims[0]*cartDims[1]) + { + mult = std::min(mult, transMult.getMultiplier(cartElemIdx, FaceDir::ZPlus)); + } + } + + trans *= mult; + applyMultipliers_(trans, outsideFaceIdx, outsideCartElemIdx, transMult); + } + 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 + auto elemIt = gridView_.template begin(); + const auto& elemEndIt = gridView_.template end(); + + for (; elemIt != elemEndIt; ++elemIt) { + const auto& elem = *elemIt; + auto isIt = gridView_.ibegin(elem); + const auto& isEndIt = gridView_.iend(elem); + for (; isIt != isEndIt; ++ isIt) { + // store intersection, this might be costly + const auto& intersection = *isIt; + if (!intersection.neighbor()) + continue; // intersection is on the domain boundary + + // In the EclState TRANX[c1] is transmissibility in X+ + // direction. Ordering of compressed (c1,c2) and cartesian index + // (gc1, gc2) is coherent (c1 < c2 <=> gc1 < gc2). This also + // holds for the global grid. While distributing changes the + // order of the local indices, the transmissibilities are still + // stored at the cell with the lower global cartesian index as + // the fieldprops are communicated by the grid. + unsigned c1 = elemMapper.index(intersection.inside()); + unsigned c2 = elemMapper.index(intersection.outside()); + int gc1 = cartMapper_.cartesianIndex(c1); + int gc2 = cartMapper_.cartesianIndex(c2); + + if (gc1 > gc2) + continue; // we only need to handle each connection once, thank you. + + auto isID = isId(c1, c2); + + if (gc2 - gc1 == 1 && cartDims[0] > 1) { + if (is_tran[0]) + // set simulator internal transmissibilities to values from inputTranx + trans[0][c1] = trans_[isID]; + } + else if (gc2 - gc1 == cartDims[0] && cartDims[1] > 1) { + if (is_tran[1]) + // set simulator internal transmissibilities to values from inputTrany + trans[1][c1] = trans_[isID]; + } + else if (gc2 - gc1 == cartDims[0]*cartDims[1]) { + if (is_tran[2]) + // set simulator internal transmissibilities to values from inputTranz + trans[2][c1] = trans_[isID]; + } + //else.. We don't support modification of NNC at the moment. + } + } + + return trans; +} + +template +void 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 + auto elemIt = gridView_.template begin(); + const auto& elemEndIt = gridView_.template end(); + + for (; elemIt != elemEndIt; ++elemIt) { + const auto& elem = *elemIt; + auto isIt = gridView_.ibegin(elem); + const auto& isEndIt = gridView_.iend(elem); + for (; isIt != isEndIt; ++ isIt) { + // store intersection, this might be costly + const auto& intersection = *isIt; + if (!intersection.neighbor()) + continue; // intersection is on the domain boundary + + // In the EclState TRANX[c1] is transmissibility in X+ + // direction. Ordering of compressed (c1,c2) and cartesian index + // (gc1, gc2) is coherent (c1 < c2 <=> gc1 < gc2). This also + // holds for the global grid. While distributing changes the + // order of the local indices, the transmissibilities are still + // stored at the cell with the lower global cartesian index as + // the fieldprops are communicated by the grid. + unsigned c1 = elemMapper.index(intersection.inside()); + unsigned c2 = elemMapper.index(intersection.outside()); + int gc1 = cartMapper_.cartesianIndex(c1); + int gc2 = cartMapper_.cartesianIndex(c2); + if (gc1 > gc2) + continue; // we only need to handle each connection once, thank you. + + auto isID = isId(c1, c2); + + if (gc2 - gc1 == 1 && cartDims[0] > 1) { + if (is_tran[0]) + // set simulator internal transmissibilities to values from inputTranx + trans_[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::vector& 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 low = cartesianToCompressed[c1]; + auto high = cartesianToCompressed[c2]; + + 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::vector& globalToLocal) +{ + const auto& nnc_input = eclState_.getInputNNC(); + const auto& editNnc = nnc_input.edit(); + if (editNnc.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 make_warning = [&format_ijk] (const KeywordLocation& location, const NNCdata& nnc) -> std::string { + return fmt::format("Problem with EDITNNC keyword\n" + "In {} line {} \n" + "No NNC defined for connection {} -> {}", location.filename, location.lineno, format_ijk(nnc.cell1), format_ijk(nnc.cell2)); + + }; + + // 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 = editNnc.begin(); + auto end = editNnc.end(); + std::size_t warning_count = 0; + while (nnc != end) { + auto c1 = nnc->cell1; + auto c2 = nnc->cell2; + auto low = globalToLocal[c1]; + auto high = globalToLocal[c2]; + if (low > high) + std::swap(low, high); + + auto candidate = trans_.find(isId(low, high)); + if (candidate == trans_.end()) { + const auto& location = nnc_input.edit_location( *nnc ); + auto warning = make_warning(location, *nnc); + OpmLog::warning("EDITNNC", warning); + ++nnc; + warning_count++; + } + else { + // NNC exists + while (nnc!= end && c1==nnc->cell1 && c2==nnc->cell2) { + candidate->second *= nnc->trans; + ++nnc; + } + } + } + + if (warning_count > 0) { + auto warning = fmt::format("Problems with EDITNNC keyword\n" + "A total of {} connections not defined in grid", 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>>, + Dune::MultipleCodimMultipleGeomTypeMapper>>, + double, + true, + true>; +template class EclTransmissibility>>, + Dune::MultipleCodimMultipleGeomTypeMapper>>, + double, + false, + true>; +#else +template class EclTransmissibility>, + Dune::MultipleCodimMultipleGeomTypeMapper>, Dune::Impl::MCMGFailLayout>, + double, + true, + true>; +template class EclTransmissibility>, + Dune::MultipleCodimMultipleGeomTypeMapper>, Dune::Impl::MCMGFailLayout>, + double, + false, + true>; +#endif + +template class EclTransmissibility, + Dune::GridView>, + Dune::MultipleCodimMultipleGeomTypeMapper>, Dune::Impl::MCMGFailLayout>, + double, + false, + true>; + +} // namespace Opm diff --git a/ebos/ecltransmissibility.hh b/ebos/ecltransmissibility.hh index 4d1cd0ecf..a443d698e 100644 --- a/ebos/ecltransmissibility.hh +++ b/ebos/ecltransmissibility.hh @@ -28,444 +28,40 @@ #ifndef EWOMS_ECL_TRANSMISSIBILITY_HH #define EWOMS_ECL_TRANSMISSIBILITY_HH -#include -#include - -#include -#include -#include -#include -#include - - -#include +#include #include -#include - -#include #include #include -#include #include +#include +#include #include #include namespace Opm { -/*! - * \ingroup EclBlackOilSimulator - * - * \brief This class calculates the transmissibilites for grid faces according to the - * Eclipse Technical Description. - */ -template -class EclTransmissibility -{ - using Grid = GetPropType; - using GridView = GetPropType; - using Scalar = GetPropType; - using Vanguard = GetPropType; - using ElementMapper = GetPropType; - using Intersection = typename GridView::Intersection; - - static constexpr bool enableEnergy = getPropValue(); - static constexpr bool enableDiffusion = getPropValue(); +class EclipseState; +struct NNCdata; +class TransMult; +template +class EclTransmissibility { // Grid and world dimension enum { dimWorld = GridView::dimensionworld }; - - typedef Dune::FieldMatrix DimMatrix; - typedef Dune::FieldVector DimVector; - - static const unsigned elemIdxShift = 32; // bits - public: - EclTransmissibility(const Vanguard& vanguard) - : vanguard_(vanguard) - { - const Opm::UnitSystem& unitSystem = vanguard_.eclState().getDeckUnitSystem(); - transmissibilityThreshold_ = unitSystem.parse("Transmissibility").getSIScaling() * 1e-6; - } + using DimMatrix = Dune::FieldMatrix; + using DimVector = Dune::FieldVector; - /*! - * \brief Actually compute the transmissibilty over a face as a pre-compute step. - * - * This code actually uses the direction specific "centroids" of - * each element. These "centroids" are _not_ the identical - * barycenter of the element, but the middle of the centers of the - * faces of the logical Cartesian cells, i.e., the centers of the - * faces of the reference elements. We do things this way because - * the barycenter of the element can be located outside of the - * element for sufficiently "ugly" (i.e., thin and "non-flat") - * elements which in turn leads to quite wrong - * permeabilities. This approach is probably not always correct - * either but at least it seems to be much better. - */ - void finishInit() - { update(true); } - - - /*! - * \brief Compute all transmissibilities - * - * \param global If true, update is called on all processes - * Also, this updates the "thermal half transmissibilities" if energy is enabled. - */ - void update(bool global) - { - const auto& gridView = vanguard_.gridView(); - const auto& cartMapper = vanguard_.cartesianIndexMapper(); - const auto& eclState = vanguard_.eclState(); - const auto& cartDims = cartMapper.cartesianDimensions(); - auto& transMult = eclState.getTransMult(); - const auto& comm = vanguard_.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() && enableDiffusion; - unsigned numElements = elemMapper.size(); - - extractPermeability_(); - - // calculate the axis specific centroids of all elements - std::array, dimWorld> axisCentroids; - - for (unsigned dimIdx = 0; dimIdx < dimWorld; ++dimIdx) - axisCentroids[dimIdx].resize(numElements); - - const std::vector& centroids = vanguard_.cellCentroids(); - - auto elemIt = gridView.template begin(); - const auto& elemEndIt = gridView.template end(); - size_t centroidIdx = 0; - for (; elemIt != elemEndIt; ++elemIt, ++centroidIdx) { - const auto& elem = *elemIt; - 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; - if (vanguard_.gridView().comm().rank() == 0) { - const auto& eclGrid = eclState.getInputGrid(); - unsigned cartesianCellIdx = cartMapper.cartesianIndex(elemIdx); - centroid = eclGrid.getCellCenter(cartesianCellIdx); - } else - std::copy(centroids.begin() + centroidIdx * dimWorld, - centroids.begin() + (centroidIdx + 1) * dimWorld, - centroid.begin()); - - 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 constexpr (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() == Opm::PinchMode::ModeEnum::ALL; - } - if (global && comm.size() > 1) { - comm.broadcast(&useSmallestMultiplier, 1, 0); - } - - // compute the transmissibilities for all intersections - elemIt = gridView.template begin(); - for (; elemIt != elemEndIt; ++elemIt) { - const auto& elem = *elemIt; - 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 constexpr (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) - Opm::FaceDir::DirEnum faceDir; - switch (insideFaceIdx) { - case 0: - case 1: - faceDir = Opm::FaceDir::XPlus; - break; - - case 2: - case 3: - faceDir = Opm::FaceDir::YPlus; - break; - - case 4: - case 5: - faceDir = Opm::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 constexpr (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 - const size_t cartesianSize = cartMapper.cartesianSize(); - // reserve memory - std::vector globalToLocal(cartesianSize, -1); - - // loop over all elements (global grid) and store Cartesian index - elemIt = vanguard_.grid().leafGridView().template begin<0>(); - - for (; elemIt != elemEndIt; ++elemIt) { - int elemIdx = elemMapper.index(*elemIt); - int cartElemIdx = vanguard_.cartesianIndexMapper().cartesianIndex(elemIdx); - globalToLocal[cartElemIdx] = elemIdx; - } - applyEditNncToGridTrans_(globalToLocal); - applyNncToGridTrans_(globalToLocal); - - //remove very small non-neighbouring transmissibilities - removeSmallNonCartesianTransmissibilities_(); - } + EclTransmissibility(const EclipseState& eclState, + const GridView& gridView, + const Dune::CartesianIndexMapper& cartMapper, + const Grid& grid, + const std::vector& centroids); /*! * \brief Return the permeability for an element. @@ -476,14 +72,12 @@ public: /*! * \brief Return the transmissibility for the intersection between two elements. */ - Scalar transmissibility(unsigned elemIdx1, unsigned elemIdx2) const - { return trans_.at(isId_(elemIdx1, elemIdx2)); } + Scalar transmissibility(unsigned elemIdx1, unsigned elemIdx2) const; /*! * \brief Return the transmissibility for a given boundary segment. */ - Scalar transmissibilityBoundary(unsigned elemIdx, unsigned boundaryFaceIdx) const - { return transBoundary_.at(std::make_pair(elemIdx, boundaryFaceIdx)); } + Scalar transmissibilityBoundary(unsigned elemIdx, unsigned boundaryFaceIdx) const; /*! * \brief Return the thermal "half transmissibility" for the intersection between two @@ -499,46 +93,44 @@ public: * is the outer unit normal, and d is the distance between the center of the inside * cell and the center of the intersection. */ - Scalar thermalHalfTrans(unsigned insideElemIdx, unsigned outsideElemIdx) const - { return thermalHalfTrans_->at(directionalIsId_(insideElemIdx, outsideElemIdx)); } + Scalar thermalHalfTrans(unsigned insideElemIdx, unsigned outsideElemIdx) const; - Scalar thermalHalfTransBoundary(unsigned insideElemIdx, unsigned boundaryFaceIdx) const - { return thermalHalfTransBoundary_.at(std::make_pair(insideElemIdx, boundaryFaceIdx)); } + Scalar thermalHalfTransBoundary(unsigned insideElemIdx, unsigned boundaryFaceIdx) const; /*! * \brief Return the diffusivity for the intersection between two elements. */ - Scalar diffusivity(unsigned elemIdx1, unsigned elemIdx2) const - { - if(diffusivity_->empty()) - return 0.0; + Scalar diffusivity(unsigned elemIdx1, unsigned elemIdx2) const; - return diffusivity_->at(isId_(elemIdx1, elemIdx2)); + /*! + * \brief Actually compute the transmissibility over a face as a pre-compute step. + * + * This code actually uses the direction specific "centroids" of + * each element. These "centroids" are _not_ the identical + * barycenter of the element, but the middle of the centers of the + * faces of the logical Cartesian cells, i.e., the centers of the + * faces of the reference elements. We do things this way because + * the barycenter of the element can be located outside of the + * element for sufficiently "ugly" (i.e., thin and "non-flat") + * elements which in turn leads to quite wrong + * permeabilities. This approach is probably not always correct + * either but at least it seems to be much better. + */ + void finishInit() + { update(true); } - } + /*! + * \brief Compute all transmissibilities + * + * \param global If true, update is called on all processes + * Also, this updates the "thermal half transmissibilities" if energy is enabled. + */ + void update(bool global); -private: +protected: + void updateFromEclState_(bool global); - void removeSmallNonCartesianTransmissibilities_() - { - const auto& cartMapper = vanguard_.cartesianIndexMapper(); - 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; - } - } - } + void removeSmallNonCartesianTransmissibilities_(); /// \brief Apply the Multipliers for the case PINCH(4)==TOPBOT /// @@ -548,221 +140,34 @@ private: unsigned outsideFaceIdx, unsigned insideCartElemIdx, unsigned outsideCartElemIdx, - const Opm::TransMult& transMult, + 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 - Scalar mult = transMult.getMultiplier(lastCartElemIdx , Opm::FaceDir::ZPlus); - - if ( !pinchTop ) - { - // pick the smallest multiplier for Z+ while looking down the pillar until reaching the other end of the connection - // While Z- is not all used here. - for(auto cartElemIdx = insideCartElemIdx; cartElemIdx < lastCartElemIdx; - cartElemIdx += cartDims[0]*cartDims[1]) - { - mult = std::min(mult, transMult.getMultiplier(cartElemIdx, Opm::FaceDir::ZPlus)); - } - } - - trans *= mult; - applyMultipliers_(trans, outsideFaceIdx, outsideCartElemIdx, transMult); - } - else - { - applyMultipliers_(trans, insideFaceIdx, insideCartElemIdx, transMult); - applyMultipliers_(trans, outsideFaceIdx, outsideCartElemIdx, transMult); - } - } - - void updateFromEclState_(bool global) - { - const FieldPropsManager* fp = - (global) ? &(vanguard_.eclState().fieldProps()) : - &(vanguard_.eclState().globalFieldProps()); - - std::array is_tran {fp->tran_active("TRANX"), - fp->tran_active("TRANY"), - fp->tran_active("TRANZ")}; - - if( !(is_tran[0] ||is_tran[1] || is_tran[2]) ) - { - // Skip unneeded expensive traversals - return; - } - - std::array keywords {"TRANX", "TRANY", "TRANZ"}; - std::array,3> trans = createTransmissibilityArrays(is_tran); - auto key = keywords.begin(); - auto perform = is_tran.begin(); - - for (auto it = trans.begin(); it != trans.end(); ++it, ++key, ++perform) - { - if(perform) - fp->apply_tran(*key, *it); - } - - resetTransmissibilityFromArrays(is_tran, trans); - } - - /// \brief overwrites calculated transmissibilities - /// - /// \param is_tran Whether TRAN{XYZ} have been modified. - /// \param trans Arrays with modified transmissibilities TRAN{XYZ} - void resetTransmissibilityFromArrays(const std::array& is_tran, - const std::array,3>& trans) - { - const auto& gridView = vanguard_.gridView(); - const auto& cartMapper = vanguard_.cartesianIndexMapper(); - const auto& cartDims = cartMapper.cartesianDimensions(); - ElementMapper elemMapper(gridView, Dune::mcmgElementLayout()); - - // compute the transmissibilities for all intersections - auto elemIt = gridView.template begin(); - const auto& elemEndIt = gridView.template end(); - - for (; elemIt != elemEndIt; ++elemIt) { - const auto& elem = *elemIt; - auto isIt = gridView.ibegin(elem); - const auto& isEndIt = gridView.iend(elem); - for (; isIt != isEndIt; ++ isIt) { - // store intersection, this might be costly - const auto& intersection = *isIt; - if (!intersection.neighbor()) - continue; // intersection is on the domain boundary - - // In the EclState TRANX[c1] is transmissibility in X+ - // direction. Ordering of compressed (c1,c2) and cartesian index - // (gc1, gc2) is coherent (c1 < c2 <=> gc1 < gc2). This also - // holds for the global grid. While distributing changes the - // order of the local indices, the transmissibilities are still - // stored at the cell with the lower global cartesian index as - // the fieldprops are communicated by the grid. - unsigned c1 = elemMapper.index(intersection.inside()); - unsigned c2 = elemMapper.index(intersection.outside()); - int gc1 = cartMapper.cartesianIndex(c1); - int gc2 = cartMapper.cartesianIndex(c2); - if (gc1 > gc2) - continue; // we only need to handle each connection once, thank you. - - auto isId = isId_(c1, c2); - - if (gc2 - gc1 == 1 && cartDims[0] > 1) { - if (is_tran[0]) - // set simulator internal transmissibilities to values from inputTranx - trans_[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. - } - } - } + bool pinchTop); /// \brief Creates TRANS{XYZ} arrays for modification by FieldProps data /// /// \param is_tran Whether TRAN{XYZ} will be modified. If entry is false the array will be empty /// \returns an array of vector (TRANX, TRANY, TRANZ} std::array,3> - createTransmissibilityArrays(const std::array& is_tran) - { - const auto& gridView = vanguard_.gridView(); - const auto& cartMapper = vanguard_.cartesianIndexMapper(); - const auto& cartDims = cartMapper.cartesianDimensions(); - ElementMapper elemMapper(gridView, Dune::mcmgElementLayout()); - - auto numElem = vanguard_.gridView().size(/*codim=*/0); - std::array,3> trans = - { std::vector(is_tran[0] ? numElem : 0, 0), - std::vector(is_tran[1] ? numElem : 0, 0), - std::vector(is_tran[2] ? numElem : 0, 0)}; - - // compute the transmissibilities for all intersections - auto elemIt = gridView.template begin(); - const auto& elemEndIt = gridView.template end(); - - for (; elemIt != elemEndIt; ++elemIt) { - const auto& elem = *elemIt; - auto isIt = gridView.ibegin(elem); - const auto& isEndIt = gridView.iend(elem); - for (; isIt != isEndIt; ++ isIt) { - // store intersection, this might be costly - const auto& intersection = *isIt; - if (!intersection.neighbor()) - continue; // intersection is on the domain boundary - - // In the EclState TRANX[c1] is transmissibility in X+ - // direction. Ordering of compressed (c1,c2) and cartesian index - // (gc1, gc2) is coherent (c1 < c2 <=> gc1 < gc2). This also - // holds for the global grid. While distributing changes the - // order of the local indices, the transmissibilities are still - // stored at the cell with the lower global cartesian index as - // the fieldprops are communicated by the grid. - unsigned c1 = elemMapper.index(intersection.inside()); - unsigned c2 = elemMapper.index(intersection.outside()); - int gc1 = cartMapper.cartesianIndex(c1); - int gc2 = cartMapper.cartesianIndex(c2); - - if (gc1 > gc2) - continue; // we only need to handle each connection once, thank you. - - auto isId = isId_(c1, c2); - - if (gc2 - gc1 == 1 && cartDims[0] > 1) { - if (is_tran[0]) - // set simulator internal transmissibilities to values from inputTranx - trans[0][c1] = trans_[isId]; - } - else if (gc2 - gc1 == cartDims[0] && cartDims[1] > 1) { - if (is_tran[1]) - // set simulator internal transmissibilities to values from inputTrany - trans[1][c1] = trans_[isId]; - } - else if (gc2 - gc1 == cartDims[0]*cartDims[1]) { - if (is_tran[2]) - // set simulator internal transmissibilities to values from inputTranz - trans[2][c1] = trans_[isId]; - } - //else.. We don't support modification of NNC at the moment. - } - } - - return trans; - } + createTransmissibilityArrays_(const std::array& is_tran); + /// \brief overwrites calculated transmissibilities + /// + /// \param is_tran Whether TRAN{XYZ} have been modified. + /// \param trans Arrays with modified transmissibilities TRAN{XYZ} + void resetTransmissibilityFromArrays_(const std::array& is_tran, + const std::array,3>& trans); template void computeFaceProperties(const Intersection& intersection, - const int insideElemIdx OPM_UNUSED, - const int insideFaceIdx OPM_UNUSED, - const int outsideElemIdx OPM_UNUSED, - const int outsideFaceIdx OPM_UNUSED, + 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(); - } + /*isCpGrid=*/std::false_type) const; template void computeFaceProperties(const Intersection& intersection, @@ -773,13 +178,7 @@ private: DimVector& faceCenterInside, DimVector& faceCenterOutside, DimVector& faceAreaNormal, - /*isCpGrid=*/std::true_type) const - { - int faceIdx = intersection.id(); - faceCenterInside = vanguard_.grid().faceCenterEcl(insideElemIdx, insideFaceIdx); - faceCenterOutside = vanguard_.grid().faceCenterEcl(outsideElemIdx, outsideFaceIdx); - faceAreaNormal = vanguard_.grid().faceAreaNormalEcl(faceIdx); - } + /*isCpGrid=*/std::true_type) const; /* * \brief Applies additional transmissibilities specified via NNC keyword. @@ -795,309 +194,51 @@ private: * and the second the NNCs not resembled by faces of the grid. NNCs specified for * inactive cells are omitted in these vectors. */ - std::tuple, std::vector > - applyNncToGridTrans_(const std::vector& cartesianToCompressed) - { - // First scale NNCs with EDITNNC. - std::vector unprocessedNnc; - std::vector processedNnc; - const auto& nnc_input = vanguard_.eclState().getInputNNC().input(); - if (nnc_input.empty()) - return make_tuple(processedNnc, unprocessedNnc); - - for (const auto& nncEntry : nnc_input) { - auto c1 = nncEntry.cell1; - auto c2 = nncEntry.cell2; - auto low = cartesianToCompressed[c1]; - auto high = cartesianToCompressed[c2]; - - 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 <<")"; - Opm::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 make_tuple(processedNnc, unprocessedNnc); - } - + std::tuple, std::vector> + applyNncToGridTrans_(const std::vector& cartesianToCompressed); /// \brief Multiplies the grid transmissibilities according to EDITNNC. - void applyEditNncToGridTrans_(const std::vector& globalToLocal) - { - const auto& nnc_input = vanguard_.eclState().getInputNNC(); - const auto& editNnc = nnc_input.edit(); - if (editNnc.empty()) - return; - const auto& cartMapper = vanguard_.cartesianIndexMapper(); - const auto& cartDims = cartMapper.cartesianDimensions(); + void applyEditNncToGridTrans_(const std::vector& globalToLocal); - 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]; + void extractPermeability_(); - return fmt::format("({},{},{})", i + 1,j + 1,k + 1); - }; - - auto make_warning = [&format_ijk] (const Opm::KeywordLocation& location, const Opm::NNCdata& nnc) -> std::string { - return fmt::format("Problem with EDITNNC keyword\n" - "In {} line {} \n" - "No NNC defined for connection {} -> {}", location.filename, location.lineno, format_ijk(nnc.cell1), format_ijk(nnc.cell2)); - - }; - - // 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 = editNnc.begin(); - auto end = editNnc.end(); - std::size_t warning_count = 0; - while (nnc != end) { - auto c1 = nnc->cell1; - auto c2 = nnc->cell2; - auto low = globalToLocal[c1]; - auto high = globalToLocal[c2]; - if (low > high) - std::swap(low, high); - - auto candidate = trans_.find(isId_(low, high)); - if (candidate == trans_.end()) { - const auto& location = nnc_input.edit_location( *nnc ); - auto warning = make_warning(location, *nnc); - Opm::OpmLog::warning("EDITNNC", warning); - ++nnc; - warning_count++; - } - else { - // NNC exists - while (nnc!= end && c1==nnc->cell1 && c2==nnc->cell2) { - candidate->second *= nnc->trans; - ++nnc; - } - } - } - - if (warning_count > 0) { - auto warning = fmt::format("Problems with EDITNNC keyword\n" - "A total of {} connections not defined in grid", warning_count); - Opm::OpmLog::warning(warning); - } - } - - void extractPermeability_() - { - unsigned numElem = vanguard_.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 = vanguard_.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)"); - } - - void 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 = vanguard_.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)"); - } - - std::uint64_t isId_(std::uint32_t elemIdx1, std::uint32_t elemIdx2) const - { - std::uint32_t elemAIdx = std::min(elemIdx1, elemIdx2); - std::uint64_t elemBIdx = std::max(elemIdx1, elemIdx2); - - return (elemBIdx< isIdReverse_(const std::uint64_t& id) const - { - // 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) const - { - return (std::uint64_t(elemIdx1)<= 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(); - } + const DimMatrix& perm) const; void 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(); - } + const Scalar& poro) const; DimVector 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; - } + const std::array, dimWorld>& axisCentroids) const; void applyMultipliers_(Scalar& trans, unsigned faceIdx, unsigned cartElemIdx, - const Opm::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, Opm::FaceDir::XMinus); - break; - case 1: // right - trans *= transMult.getMultiplier(cartElemIdx, Opm::FaceDir::XPlus); - break; - - case 2: // front - trans *= transMult.getMultiplier(cartElemIdx, Opm::FaceDir::YMinus); - break; - case 3: // back - trans *= transMult.getMultiplier(cartElemIdx, Opm::FaceDir::YPlus); - break; - - case 4: // bottom - trans *= transMult.getMultiplier(cartElemIdx, Opm::FaceDir::ZMinus); - break; - case 5: // top - trans *= transMult.getMultiplier(cartElemIdx, Opm::FaceDir::ZPlus); - break; - } - } + const TransMult& transMult) const; void 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; + const std::vector& ntg) const; - case 2: // front - trans *= ntg[elemIdx]; - break; - case 3: // back - trans *= ntg[elemIdx]; - break; - - // NTG does not apply to top and bottom faces - } - } - - const Vanguard& vanguard_; - Scalar transmissibilityThreshold_; std::vector permeability_; std::vector porosity_; std::unordered_map trans_; + const EclipseState& eclState_; + const GridView& gridView_; + const Dune::CartesianIndexMapper& cartMapper_; + const Grid& grid_; + const std::vector& centroids_; + Scalar transmissibilityThreshold_; std::map, Scalar> transBoundary_; std::map, Scalar> thermalHalfTransBoundary_; Opm::ConditionalStorage ElementMapper; ElementMapper globalElemMapper(globalGridView, Dune::mcmgElementLayout()); - const EclTransmissibility* globalTrans; + using TransmissibilityType = typename Vanguard::TransmissibilityType; + + const TransmissibilityType* globalTrans; if (!collectToIORank_.isParallel()) { @@ -622,7 +624,8 @@ private: typedef Dune::MultipleCodimMultipleGeomTypeMapper ElementMapper; ElementMapper globalElemMapper(globalGridView, Dune::mcmgElementLayout()); - const EclTransmissibility* globalTrans; + using TransmissibilityType = typename Vanguard::TransmissibilityType; + const TransmissibilityType* globalTrans; if (!collectToIORank_.isParallel()) { // in the sequential case we must use the transmissibilites defined by // the problem. (because in the sequential case, the grid manager does