23#ifndef OPM_TRANSMISSIBILITY_IMPL_HPP
24#define OPM_TRANSMISSIBILITY_IMPL_HPP
26#include <dune/common/version.hh>
27#include <dune/grid/common/mcmgmapper.hh>
29#include <opm/grid/CpGrid.hpp>
31#include <opm/common/OpmLog/KeywordLocation.hpp>
32#include <opm/input/eclipse/EclipseState/EclipseState.hpp>
33#include <opm/input/eclipse/EclipseState/Grid/FaceDir.hpp>
34#include <opm/input/eclipse/EclipseState/Grid/FieldPropsManager.hpp>
35#include <opm/input/eclipse/EclipseState/Grid/TransMult.hpp>
36#include <opm/input/eclipse/Units/Units.hpp>
40#include <fmt/format.h>
49#include <initializer_list>
62 std::uint64_t
isId(std::uint32_t elemIdx1, std::uint32_t elemIdx2)
64 std::uint32_t elemAIdx = std::min(elemIdx1, elemIdx2);
65 std::uint64_t elemBIdx = std::max(elemIdx1, elemIdx2);
70 std::pair<std::uint32_t, std::uint32_t>
isIdReverse(
const std::uint64_t&
id)
75 std::uint32_t elemAIdx = id;
76 std::uint32_t elemBIdx = (
id - elemAIdx) >>
elemIdxShift;
78 return std::make_pair(elemAIdx, elemBIdx);
83 return (std::uint64_t(elemIdx1)<<
elemIdxShift) + elemIdx2;
87template<
class Gr
id,
class Gr
idView,
class ElementMapper,
class CartesianIndexMapper,
class Scalar>
90 const GridView& gridView,
91 const CartesianIndexMapper& cartMapper,
93 std::function<std::array<double,dimWorld>(
int)> centroids,
95 bool enableDiffusivity,
96 bool enableDispersivity)
99 , cartMapper_(cartMapper)
101 , centroids_(centroids)
102 , enableEnergy_(enableEnergy)
103 , enableDiffusivity_(enableDiffusivity)
104 , enableDispersivity_(enableDispersivity)
105 , lookUpData_(gridView)
106 , lookUpCartesianData_(gridView, cartMapper)
108 const UnitSystem& unitSystem =
eclState_.getDeckUnitSystem();
112template<
class Gr
id,
class Gr
idView,
class ElementMapper,
class CartesianIndexMapper,
class Scalar>
119template<
class Gr
id,
class Gr
idView,
class ElementMapper,
class CartesianIndexMapper,
class Scalar>
123 return transBoundary_.at(std::make_pair(elemIdx, boundaryFaceIdx));
126template<
class Gr
id,
class Gr
idView,
class ElementMapper,
class CartesianIndexMapper,
class Scalar>
133template<
class Gr
id,
class Gr
idView,
class ElementMapper,
class CartesianIndexMapper,
class Scalar>
137 return thermalHalfTransBoundary_.at(std::make_pair(insideElemIdx, boundaryFaceIdx));
140template<
class Gr
id,
class Gr
idView,
class ElementMapper,
class CartesianIndexMapper,
class Scalar>
142diffusivity(
unsigned elemIdx1,
unsigned elemIdx2)
const
144 if (diffusivity_.empty())
150template<
class Gr
id,
class Gr
idView,
class ElementMapper,
class CartesianIndexMapper,
class Scalar>
154 if (dispersivity_.empty())
160template<
class Gr
id,
class Gr
idView,
class ElementMapper,
class CartesianIndexMapper,
class Scalar>
162update(
bool global,
const std::function<
unsigned int(
unsigned int)>& map,
const bool applyNncMultregT)
164 const auto& cartDims = cartMapper_.cartesianDimensions();
165 const auto& transMult = eclState_.getTransMult();
166 const auto& comm = gridView_.comm();
167 ElementMapper elemMapper(gridView_, Dune::mcmgElementLayout());
169 unsigned numElements = elemMapper.size();
171 const std::vector<double>& ntg = this->lookUpData_.assignFieldPropsDoubleOnLeaf(eclState_.fieldProps(),
"NTG");
172 const bool updateDiffusivity = eclState_.getSimulationConfig().isDiffusive();
173 const bool updateDispersivity = eclState_.getSimulationConfig().rock_config().dispersion();
175 const bool disableNNC = eclState_.getSimulationConfig().useNONNC();
178 extractPermeability_(map);
180 extractPermeability_();
183 std::array<std::vector<DimVector>, dimWorld> axisCentroids;
185 for (
unsigned dimIdx = 0; dimIdx < dimWorld; ++dimIdx)
186 axisCentroids[dimIdx].resize(numElements);
188 for (
const auto& elem : elements(gridView_)) {
189 unsigned elemIdx = elemMapper.index(elem);
194 std::array<double, dimWorld> centroid = centroids_(elemIdx);
196 for (
unsigned axisIdx = 0; axisIdx < dimWorld; ++axisIdx)
197 for (
unsigned dimIdx = 0; dimIdx < dimWorld; ++dimIdx)
198 axisCentroids[axisIdx][elemIdx][dimIdx] = centroid[dimIdx];
206 trans_.reserve(numElements*3*1.05);
208 transBoundary_.clear();
212 thermalHalfTrans_.clear();
213 thermalHalfTrans_.reserve(numElements*6*1.05);
215 thermalHalfTransBoundary_.clear();
219 if (updateDiffusivity) {
220 diffusivity_.clear();
221 diffusivity_.reserve(numElements*3*1.05);
226 if (updateDispersivity) {
227 dispersivity_.clear();
228 dispersivity_.reserve(numElements*3*1.05);
229 extractDispersion_();
235 bool useSmallestMultiplier;
237 if (comm.rank() == 0) {
238 const auto& eclGrid = eclState_.getInputGrid();
239 pinchActive = eclGrid.isPinchActive();
240 useSmallestMultiplier = eclGrid.getMultzOption() == PinchMode::ALL;
242 if (global && comm.size() > 1) {
243 comm.broadcast(&useSmallestMultiplier, 1, 0);
244 comm.broadcast(&pinchActive, 1, 0);
248 for (
const auto& elem : elements(gridView_)) {
249 unsigned elemIdx = elemMapper.index(elem);
251 auto isIt = gridView_.ibegin(elem);
252 const auto& isEndIt = gridView_.iend(elem);
253 unsigned boundaryIsIdx = 0;
254 for (; isIt != isEndIt; ++ isIt) {
256 const auto& intersection = *isIt;
259 if (intersection.boundary()) {
261 const auto& geometry = intersection.geometry();
262 const auto& faceCenterInside = geometry.center();
264 auto faceAreaNormal = intersection.centerUnitOuterNormal();
265 faceAreaNormal *= geometry.volume();
267 Scalar transBoundaryIs;
268 computeHalfTrans_(transBoundaryIs,
270 intersection.indexInInside(),
271 distanceVector_(faceCenterInside,
272 intersection.indexInInside(),
275 permeability_[elemIdx]);
280 unsigned insideCartElemIdx = cartMapper_.cartesianIndex(elemIdx);
281 applyMultipliers_(transBoundaryIs, intersection.indexInInside(), insideCartElemIdx, transMult);
282 transBoundary_[std::make_pair(elemIdx, boundaryIsIdx)] = transBoundaryIs;
287 Scalar transBoundaryEnergyIs;
288 computeHalfDiffusivity_(transBoundaryEnergyIs,
290 distanceVector_(faceCenterInside,
291 intersection.indexInInside(),
295 thermalHalfTransBoundary_[std::make_pair(elemIdx, boundaryIsIdx)] =
296 transBoundaryEnergyIs;
303 if (!intersection.neighbor()) {
310 const auto& outsideElem = intersection.outside();
311 unsigned outsideElemIdx = elemMapper.index(outsideElem);
315 unsigned insideCartElemIdx = this-> lookUpCartesianData_.template getFieldPropCartesianIdx<Grid>(elemIdx);
316 unsigned outsideCartElemIdx = this-> lookUpCartesianData_.template getFieldPropCartesianIdx<Grid>(outsideElemIdx);
324 if (std::tie(insideCartElemIdx, elemIdx) > std::tie(outsideCartElemIdx, outsideElemIdx))
329 int insideFaceIdx = intersection.indexInInside();
330 int outsideFaceIdx = intersection.indexInOutside();
332 if (insideFaceIdx == -1) {
335 assert(outsideFaceIdx == -1);
342 if (updateDiffusivity) {
345 if (updateDispersivity) {
355 typename std::is_same<Grid, Dune::CpGrid>::type isCpGrid;
356 computeFaceProperties(intersection,
369 computeHalfTrans_(halfTrans1,
372 distanceVector_(faceCenterInside,
373 intersection.indexInInside(),
376 permeability_[elemIdx]);
377 computeHalfTrans_(halfTrans2,
380 distanceVector_(faceCenterOutside,
381 intersection.indexInOutside(),
384 permeability_[outsideElemIdx]);
386 applyNtg_(halfTrans1, insideFaceIdx, elemIdx, ntg);
387 applyNtg_(halfTrans2, outsideFaceIdx, outsideElemIdx, ntg);
392 if (std::abs(halfTrans1) < 1e-30 || std::abs(halfTrans2) < 1e-30)
396 trans = 1.0 / (1.0/halfTrans1 + 1.0/halfTrans2);
401 if (insideFaceIdx > 3){
402 auto find_layer = [&cartDims](std::size_t cell){
404 auto k = cell / cartDims[1];
407 int kup = find_layer(insideCartElemIdx);
408 int kdown=find_layer(outsideCartElemIdx);
411 assert((kup != kdown) || (insideCartElemIdx == outsideCartElemIdx));
412 if(std::abs(kup -kdown) > 1){
418 if (useSmallestMultiplier)
422 applyAllZMultipliers_(trans, insideFaceIdx, outsideFaceIdx, insideCartElemIdx,
423 outsideCartElemIdx, transMult, cartDims,
428 applyMultipliers_(trans, insideFaceIdx, insideCartElemIdx, transMult);
430 applyMultipliers_(trans, outsideFaceIdx, outsideCartElemIdx, transMult);
434 FaceDir::DirEnum faceDir;
435 switch (insideFaceIdx) {
438 faceDir = FaceDir::XPlus;
443 faceDir = FaceDir::YPlus;
448 faceDir = FaceDir::ZPlus;
452 throw std::logic_error(
"Could not determine a face direction");
455 trans *= transMult.getRegionMultiplier(insideCartElemIdx,
464 Scalar halfDiffusivity1;
465 Scalar halfDiffusivity2;
467 computeHalfDiffusivity_(halfDiffusivity1,
469 distanceVector_(faceCenterInside,
470 intersection.indexInInside(),
474 computeHalfDiffusivity_(halfDiffusivity2,
476 distanceVector_(faceCenterOutside,
477 intersection.indexInOutside(),
487 if (updateDiffusivity) {
489 Scalar halfDiffusivity1;
490 Scalar halfDiffusivity2;
492 computeHalfDiffusivity_(halfDiffusivity1,
494 distanceVector_(faceCenterInside,
495 intersection.indexInInside(),
499 computeHalfDiffusivity_(halfDiffusivity2,
501 distanceVector_(faceCenterOutside,
502 intersection.indexInOutside(),
505 porosity_[outsideElemIdx]);
507 applyNtg_(halfDiffusivity1, insideFaceIdx, elemIdx, ntg);
508 applyNtg_(halfDiffusivity2, outsideFaceIdx, outsideElemIdx, ntg);
512 if (std::abs(halfDiffusivity1) < 1e-30 || std::abs(halfDiffusivity2) < 1e-30)
516 diffusivity = 1.0 / (1.0/halfDiffusivity1 + 1.0/halfDiffusivity2);
519 diffusivity_[
details::isId(elemIdx, outsideElemIdx)] = diffusivity;
523 if (updateDispersivity) {
525 Scalar halfDispersivity1;
526 Scalar halfDispersivity2;
528 computeHalfDiffusivity_(halfDispersivity1,
530 distanceVector_(faceCenterInside,
531 intersection.indexInInside(),
534 dispersion_[elemIdx]);
535 computeHalfDiffusivity_(halfDispersivity2,
537 distanceVector_(faceCenterOutside,
538 intersection.indexInOutside(),
541 dispersion_[outsideElemIdx]);
543 applyNtg_(halfDispersivity1, insideFaceIdx, elemIdx, ntg);
544 applyNtg_(halfDispersivity2, outsideFaceIdx, outsideElemIdx, ntg);
548 if (std::abs(halfDispersivity1) < 1e-30 || std::abs(halfDispersivity2) < 1e-30)
552 dispersivity = 1.0 / (1.0/halfDispersivity1 + 1.0/halfDispersivity2);
555 dispersivity_[
details::isId(elemIdx, outsideElemIdx)] = dispersivity;
561 this->updateFromEclState_(global);
564 std::unordered_map<std::size_t,int> globalToLocal;
567 for (
const auto& elem : elements(grid_.leafGridView())) {
568 int elemIdx = elemMapper.index(elem);
569 int cartElemIdx = cartMapper_.cartesianIndex(elemIdx);
570 globalToLocal[cartElemIdx] = elemIdx;
574 this->applyEditNncToGridTrans_(globalToLocal);
575 this->applyNncToGridTrans_(globalToLocal);
576 this->applyEditNncrToGridTrans_(globalToLocal);
577 if (applyNncMultregT) {
578 this->applyNncMultreg_(globalToLocal);
584 this->removeNonCartesianTransmissibilities_(disableNNC);
587template<
class Gr
id,
class Gr
idView,
class ElementMapper,
class CartesianIndexMapper,
class Scalar>
591 unsigned numElem = gridView_.size(0);
592 permeability_.resize(numElem);
598 const auto& fp = eclState_.fieldProps();
599 if (fp.has_double(
"PERMX")) {
600 const std::vector<double>& permxData = this-> lookUpData_.assignFieldPropsDoubleOnLeaf(fp,
"PERMX");
602 std::vector<double> permyData;
603 if (fp.has_double(
"PERMY"))
604 permyData = this-> lookUpData_.assignFieldPropsDoubleOnLeaf(fp,
"PERMY");
606 permyData = permxData;
608 std::vector<double> permzData;
609 if (fp.has_double(
"PERMZ"))
610 permzData = this-> lookUpData_.assignFieldPropsDoubleOnLeaf(fp,
"PERMZ");
612 permzData = permxData;
614 for (std::size_t dofIdx = 0; dofIdx < numElem; ++ dofIdx) {
615 permeability_[dofIdx] = 0.0;
616 permeability_[dofIdx][0][0] = permxData[dofIdx];
617 permeability_[dofIdx][1][1] = permyData[dofIdx];
618 permeability_[dofIdx][2][2] = permzData[dofIdx];
625 throw std::logic_error(
"Can't read the intrinsic permeability from the ecl state. "
626 "(The PERM{X,Y,Z} keywords are missing)");
629template<
class Gr
id,
class Gr
idView,
class ElementMapper,
class CartesianIndexMapper,
class Scalar>
633 unsigned numElem = gridView_.size(0);
634 permeability_.resize(numElem);
640 const auto& fp = eclState_.fieldProps();
641 if (fp.has_double(
"PERMX")) {
642 const std::vector<double>& permxData = this-> lookUpData_.assignFieldPropsDoubleOnLeaf(fp,
"PERMX");
644 std::vector<double> permyData;
645 if (fp.has_double(
"PERMY"))
646 permyData = this-> lookUpData_.assignFieldPropsDoubleOnLeaf(fp,
"PERMY");
648 permyData = permxData;
650 std::vector<double> permzData;
651 if (fp.has_double(
"PERMZ"))
652 permzData = this-> lookUpData_.assignFieldPropsDoubleOnLeaf(fp,
"PERMZ");
654 permzData = permxData;
656 for (std::size_t dofIdx = 0; dofIdx < numElem; ++ dofIdx) {
657 permeability_[dofIdx] = 0.0;
658 std::size_t inputDofIdx = map(dofIdx);
659 permeability_[dofIdx][0][0] = permxData[inputDofIdx];
660 permeability_[dofIdx][1][1] = permyData[inputDofIdx];
661 permeability_[dofIdx][2][2] = permzData[inputDofIdx];
668 throw std::logic_error(
"Can't read the intrinsic permeability from the ecl state. "
669 "(The PERM{X,Y,Z} keywords are missing)");
672template<
class Gr
id,
class Gr
idView,
class ElementMapper,
class CartesianIndexMapper,
class Scalar>
680 const auto& fp = eclState_.fieldProps();
681 if (fp.has_double(
"PORO")) {
682 porosity_ = this-> lookUpData_.assignFieldPropsDoubleOnLeaf(fp,
"PORO");
685 throw std::logic_error(
"Can't read the porosityfrom the ecl state. "
686 "(The PORO keywords are missing)");
689template<
class Gr
id,
class Gr
idView,
class ElementMapper,
class CartesianIndexMapper,
class Scalar>
693 if (!enableDispersivity_) {
694 throw std::runtime_error(
"Dispersion disabled at compile time, but the deck "
695 "contains the DISPERC keyword.");
697 const auto& fp = eclState_.fieldProps();
698 dispersion_ = this-> lookUpData_.assignFieldPropsDoubleOnLeaf(fp,
"DISPERC");
701template<
class Gr
id,
class Gr
idView,
class ElementMapper,
class CartesianIndexMapper,
class Scalar>
705 const auto& cartDims = cartMapper_.cartesianDimensions();
706 for (
auto&& trans: trans_) {
708 if (removeAll or trans.second < transmissibilityThreshold_) {
709 const auto&
id = trans.first;
711 int gc1 = std::min(cartMapper_.cartesianIndex(elements.first), cartMapper_.cartesianIndex(elements.second));
712 int gc2 = std::max(cartMapper_.cartesianIndex(elements.first), cartMapper_.cartesianIndex(elements.second));
717 if (gc2 - gc1 == 1 || gc2 - gc1 == cartDims[0] || gc2 - gc1 == cartDims[0]*cartDims[1] || gc2 - gc1 == 0)
725template<
class Gr
id,
class Gr
idView,
class ElementMapper,
class CartesianIndexMapper,
class Scalar>
728 unsigned insideFaceIdx,
729 unsigned outsideFaceIdx,
730 unsigned insideCartElemIdx,
731 unsigned outsideCartElemIdx,
732 const TransMult& transMult,
733 const std::array<int, dimWorld>& cartDims,
736 if(grid_.maxLevel()> 0) {
737 OPM_THROW(std::invalid_argument,
"MULTZ not support with LGRS, yet.");
739 if (insideFaceIdx > 3) {
740 assert(insideFaceIdx==5);
742 assert(outsideCartElemIdx >= insideCartElemIdx);
743 unsigned lastCartElemIdx;
744 if (outsideCartElemIdx == insideCartElemIdx) {
745 lastCartElemIdx = outsideCartElemIdx;
748 lastCartElemIdx = outsideCartElemIdx - cartDims[0]*cartDims[1];
751 Scalar mult = transMult.getMultiplier(lastCartElemIdx , FaceDir::ZPlus) *
752 transMult.getMultiplier(outsideCartElemIdx , FaceDir::ZMinus);
758 for(
auto cartElemIdx = insideCartElemIdx; cartElemIdx < lastCartElemIdx;)
760 auto multiplier = transMult.getMultiplier(cartElemIdx, FaceDir::ZPlus);
761 cartElemIdx += cartDims[0]*cartDims[1];
762 multiplier *= transMult.getMultiplier(cartElemIdx, FaceDir::ZMinus);
763 mult = std::min(mult, multiplier);
771 applyMultipliers_(trans, insideFaceIdx, insideCartElemIdx, transMult);
772 applyMultipliers_(trans, outsideFaceIdx, outsideCartElemIdx, transMult);
776template<
class Gr
id,
class Gr
idView,
class ElementMapper,
class CartesianIndexMapper,
class Scalar>
780 const FieldPropsManager* fp =
781 (global) ? &(eclState_.fieldProps()) :
782 &(eclState_.globalFieldProps());
784 std::array<bool,3> is_tran {fp->tran_active(
"TRANX"),
785 fp->tran_active(
"TRANY"),
786 fp->tran_active(
"TRANZ")};
788 if( !(is_tran[0] ||is_tran[1] || is_tran[2]) )
794 std::array<std::string, 3> keywords {
"TRANX",
"TRANY",
"TRANZ"};
795 std::array<std::vector<double>,3> trans = createTransmissibilityArrays_(is_tran);
796 auto key = keywords.begin();
797 auto perform = is_tran.begin();
799 for (
auto it = trans.begin(); it != trans.end(); ++it, ++key, ++perform)
802 if(grid_.maxLevel()>0) {
803 OPM_THROW(std::invalid_argument,
"Calculations on TRANX/TRANY/TRANZ arrays are not support with LGRS, yet.");
805 fp->apply_tran(*key, *it);
809 resetTransmissibilityFromArrays_(is_tran, trans);
812template<
class Gr
id,
class Gr
idView,
class ElementMapper,
class CartesianIndexMapper,
class Scalar>
813std::array<std::vector<double>,3>
817 const auto& cartDims = cartMapper_.cartesianDimensions();
818 ElementMapper elemMapper(gridView_, Dune::mcmgElementLayout());
820 auto numElem = gridView_.size(0);
821 std::array<std::vector<double>,3> trans =
822 { std::vector<double>(is_tran[0] ? numElem : 0, 0),
823 std::vector<double>(is_tran[1] ? numElem : 0, 0),
824 std::vector<double>(is_tran[2] ? numElem : 0, 0)};
827 for (
const auto& elem : elements(gridView_)) {
828 for (
const auto& intersection : intersections(gridView_, elem)) {
830 if (!intersection.neighbor())
840 unsigned c1 = elemMapper.index(intersection.inside());
841 unsigned c2 = elemMapper.index(intersection.outside());
842 int gc1 = cartMapper_.cartesianIndex(c1);
843 int gc2 = cartMapper_.cartesianIndex(c2);
854 if ((gc2 - gc1 == 1 || (gc2 == gc1 && (intersection.indexInInside() == 0 || intersection.indexInInside() == 1)))
855 && cartDims[0] > 1) {
858 trans[0][c1] = trans_[isID];
860 else if ((gc2 - gc1 == cartDims[0] || (gc2 == gc1 && (intersection.indexInInside() == 2 || intersection.indexInInside() == 3)))
861 && cartDims[1] > 1) {
864 trans[1][c1] = trans_[isID];
866 else if (gc2 - gc1 == cartDims[0]*cartDims[1] ||
867 (gc2 == gc1 && (intersection.indexInInside() == 4 || intersection.indexInInside() == 5))) {
870 trans[2][c1] = trans_[isID];
880template<
class Gr
id,
class Gr
idView,
class ElementMapper,
class CartesianIndexMapper,
class Scalar>
883 const std::array<std::vector<double>,3>& trans)
885 const auto& cartDims = cartMapper_.cartesianDimensions();
886 ElementMapper elemMapper(gridView_, Dune::mcmgElementLayout());
889 for (
const auto& elem : elements(gridView_)) {
890 for (
const auto& intersection : intersections(gridView_, elem)) {
891 if (!intersection.neighbor())
904 unsigned c1 = elemMapper.index(intersection.inside());
905 unsigned c2 = elemMapper.index(intersection.outside());
906 int gc1 = cartMapper_.cartesianIndex(c1);
907 int gc2 = cartMapper_.cartesianIndex(c2);
917 if ((gc2 - gc1 == 1 || (gc2 == gc1 && (intersection.indexInInside() == 0 || intersection.indexInInside() == 1)))
918 && cartDims[0] > 1) {
921 trans_[isID] = trans[0][c1];
923 else if ((gc2 - gc1 == cartDims[0] || (gc2 == gc1 && (intersection.indexInInside() == 2|| intersection.indexInInside() == 3)))
924 && cartDims[1] > 1) {
927 trans_[isID] = trans[1][c1];
929 else if (gc2 - gc1 == cartDims[0]*cartDims[1] ||
930 (gc2 == gc1 && (intersection.indexInInside() == 4 || intersection.indexInInside() == 5))) {
933 trans_[isID] = trans[2][c1];
941template<
class Gr
id,
class Gr
idView,
class ElementMapper,
class CartesianIndexMapper,
class Scalar>
942template<
class Intersection>
952 std::false_type)
const
955 const auto& geometry = intersection.geometry();
956 faceCenterInside = geometry.center();
957 faceCenterOutside = faceCenterInside;
959 faceAreaNormal = intersection.centerUnitOuterNormal();
960 faceAreaNormal *= geometry.volume();
963template<
class Gr
id,
class Gr
idView,
class ElementMapper,
class CartesianIndexMapper,
class Scalar>
964template<
class Intersection>
967 const int insideElemIdx,
968 const int insideFaceIdx,
969 const int outsideElemIdx,
970 const int outsideFaceIdx,
974 std::true_type)
const
976 int faceIdx = intersection.id();
978 if(grid_.maxLevel() == 0) {
979 faceCenterInside = grid_.faceCenterEcl(insideElemIdx, insideFaceIdx, intersection);
980 faceCenterOutside = grid_.faceCenterEcl(outsideElemIdx, outsideFaceIdx, intersection);
981 faceAreaNormal = grid_.faceAreaNormalEcl(faceIdx);
984 if ((intersection.inside().level() != intersection.outside().level())) {
992 const auto& parentIntersection = grid_.getParentIntersectionFromLgrBoundaryFace(intersection);
993 const auto& parentIntersectionGeometry = parentIntersection.geometry();
997 faceCenterInside = (intersection.inside().level() == 0) ? parentIntersectionGeometry.center() :
998 grid_.faceCenterEcl(insideElemIdx, insideFaceIdx, intersection);
999 faceCenterOutside = (intersection.outside().level() == 0) ? parentIntersectionGeometry.center() :
1000 grid_.faceCenterEcl(outsideElemIdx, outsideFaceIdx, intersection);
1008 faceAreaNormal = intersection.centerUnitOuterNormal();
1009 faceAreaNormal *= intersection.geometry().volume();
1012 assert(intersection.inside().level() == intersection.outside().level());
1014 faceCenterInside = grid_.faceCenterEcl(insideElemIdx, insideFaceIdx, intersection);
1015 faceCenterOutside = grid_.faceCenterEcl(outsideElemIdx, outsideFaceIdx, intersection);
1018 if (intersection.inside().level() > 0) {
1019 faceAreaNormal = intersection.centerUnitOuterNormal();
1020 faceAreaNormal *= intersection.geometry().volume();
1023 faceAreaNormal = grid_.faceAreaNormalEcl(faceIdx);
1029template<
class Gr
id,
class Gr
idView,
class ElementMapper,
class CartesianIndexMapper,
class Scalar>
1035 const auto& nnc_input = eclState_.getInputNNC().input();
1037 for (
const auto& nncEntry : nnc_input) {
1038 auto c1 = nncEntry.cell1;
1039 auto c2 = nncEntry.cell2;
1040 auto lowIt = cartesianToCompressed.find(c1);
1041 auto highIt = cartesianToCompressed.find(c2);
1042 int low = (lowIt == cartesianToCompressed.end())? -1 : lowIt->second;
1043 int high = (highIt == cartesianToCompressed.end())? -1 : highIt->second;
1046 std::swap(low, high);
1048 if (low == -1 && high == -1)
1052 if (low == -1 || high == -1) {
1054 std::ostringstream sstr;
1055 sstr <<
"NNC between active and inactive cells ("
1056 << low <<
" -> " << high <<
") with globalcell is (" << c1 <<
"->" << c2 <<
")";
1057 OpmLog::warning(sstr.str());
1063 if (candidate != trans_.end()) {
1067 candidate->second += nncEntry.trans;
1098template<
class Gr
id,
class Gr
idView,
class ElementMapper,
class CartesianIndexMapper,
class Scalar>
1102 const auto& input = eclState_.getInputNNC();
1103 applyEditNncToGridTransHelper_(globalToLocal,
"EDITNNC",
1105 [&input](
const NNCdata& nnc){
1106 return input.edit_location(nnc);},
1108 [](
double& trans,
const double& rhs){ trans *= rhs;});
1111template<
class Gr
id,
class Gr
idView,
class ElementMapper,
class CartesianIndexMapper,
class Scalar>
1115 const auto& input = eclState_.getInputNNC();
1116 applyEditNncToGridTransHelper_(globalToLocal,
"EDITNNCR",
1118 [&input](
const NNCdata& nnc){
1119 return input.editr_location(nnc);},
1121 [](
double& trans,
const double& rhs){ trans = rhs;});
1124template<
class Gr
id,
class Gr
idView,
class ElementMapper,
class CartesianIndexMapper,
class Scalar>
1127 const std::string& keyword,
1128 const std::vector<NNCdata>& nncs,
1129 const std::function<KeywordLocation(
const NNCdata&)>& getLocation,
1130 const std::function<
void(
double&,
const double&)>& apply)
1134 const auto& cartDims = cartMapper_.cartesianDimensions();
1136 auto format_ijk = [&cartDims](std::size_t cell) -> std::string {
1137 auto i = cell % cartDims[0]; cell /= cartDims[0];
1138 auto j = cell % cartDims[1];
1139 auto k = cell / cartDims[1];
1141 return fmt::format(
"({},{},{})", i + 1,j + 1,k + 1);
1144 auto print_warning = [&format_ijk, &getLocation, &keyword] (
const NNCdata& nnc) {
1145 const auto& location = getLocation( nnc );
1146 auto warning = fmt::format(
"Problem with {} keyword\n"
1148 "No NNC defined for connection {} -> {}", keyword, location.filename,
1149 location.lineno, format_ijk(nnc.cell1), format_ijk(nnc.cell2));
1150 OpmLog::warning(keyword, warning);
1156 auto nnc = nncs.begin();
1157 auto end = nncs.end();
1158 std::size_t warning_count = 0;
1159 while (nnc != end) {
1160 auto c1 = nnc->cell1;
1161 auto c2 = nnc->cell2;
1162 auto lowIt = globalToLocal.find(c1);
1163 auto highIt = globalToLocal.find(c2);
1165 if (lowIt == globalToLocal.end() || highIt == globalToLocal.end()) {
1166 print_warning(*nnc);
1172 auto low = lowIt->second, high = highIt->second;
1175 std::swap(low, high);
1178 if (candidate == trans_.end()) {
1179 print_warning(*nnc);
1185 while (nnc!= end && c1==nnc->cell1 && c2==nnc->cell2) {
1186 apply(candidate->second, nnc->trans);
1192 if (warning_count > 0) {
1193 auto warning = fmt::format(
"Problems with {} keyword\n"
1194 "A total of {} connections not defined in grid", keyword, warning_count);
1195 OpmLog::warning(warning);
1199template<
class Gr
id,
class Gr
idView,
class ElementMapper,
class CartesianIndexMapper,
class Scalar>
1202applyNncMultreg_(
const std::unordered_map<std::size_t,int>& cartesianToCompressed)
1204 const auto& inputNNC = this->eclState_.getInputNNC();
1205 const auto& transMult = this->eclState_.getTransMult();
1207 auto compressedIdx = [&cartesianToCompressed](
const std::size_t globIdx)
1209 auto ixPos = cartesianToCompressed.find(globIdx);
1210 return (ixPos == cartesianToCompressed.end()) ? -1 : ixPos->second;
1224 for (
const auto& nncList : { &NNC::input, &NNC::editr }) {
1225 for (
const auto& nncEntry : (inputNNC.*nncList)()) {
1226 const auto c1 = nncEntry.cell1;
1227 const auto c2 = nncEntry.cell2;
1229 auto low = compressedIdx(c1);
1230 auto high = compressedIdx(c2);
1232 if ((low == -1) || (high == -1)) {
1237 std::swap(low, high);
1240 auto candidate = this->trans_.find(
details::isId(low, high));
1241 if (candidate != this->trans_.end()) {
1242 candidate->second *= transMult.getRegionMultiplierNNC(c1, c2);
1248template<
class Gr
id,
class Gr
idView,
class ElementMapper,
class CartesianIndexMapper,
class Scalar>
1256 assert(faceIdx >= 0);
1257 unsigned dimIdx = faceIdx/2;
1258 assert(dimIdx < dimWorld);
1259 halfTrans = perm[dimIdx][dimIdx];
1262 for (
unsigned i = 0; i < areaNormal.size(); ++i)
1263 val += areaNormal[i]*distance[i];
1265 halfTrans *= std::abs(val);
1266 halfTrans /= distance.two_norm2();
1269template<
class Gr
id,
class Gr
idView,
class ElementMapper,
class CartesianIndexMapper,
class Scalar>
1274 const Scalar& poro)
const
1278 for (
unsigned i = 0; i < areaNormal.size(); ++i)
1279 val += areaNormal[i]*distance[i];
1281 halfDiff *= std::abs(val);
1282 halfDiff /= distance.two_norm2();
1285template<
class Gr
id,
class Gr
idView,
class ElementMapper,
class CartesianIndexMapper,
class Scalar>
1291 const std::array<std::vector<DimVector>, dimWorld>& axisCentroids)
const
1293 assert(faceIdx >= 0);
1294 unsigned dimIdx = faceIdx/2;
1295 assert(dimIdx < dimWorld);
1297 x -= axisCentroids[dimIdx][elemIdx];
1302template<
class Gr
id,
class Gr
idView,
class ElementMapper,
class CartesianIndexMapper,
class Scalar>
1306 unsigned cartElemIdx,
1307 const TransMult& transMult)
const
1314 trans *= transMult.getMultiplier(cartElemIdx, FaceDir::XMinus);
1317 trans *= transMult.getMultiplier(cartElemIdx, FaceDir::XPlus);
1321 trans *= transMult.getMultiplier(cartElemIdx, FaceDir::YMinus);
1324 trans *= transMult.getMultiplier(cartElemIdx, FaceDir::YPlus);
1328 trans *= transMult.getMultiplier(cartElemIdx, FaceDir::ZMinus);
1331 trans *= transMult.getMultiplier(cartElemIdx, FaceDir::ZPlus);
1336template<
class Gr
id,
class Gr
idView,
class ElementMapper,
class CartesianIndexMapper,
class Scalar>
1341 const std::vector<double>& ntg)
const
1348 trans *= ntg[elemIdx];
1351 trans *= ntg[elemIdx];
1355 trans *= ntg[elemIdx];
1358 trans *= ntg[elemIdx];
void applyMultipliers_(Scalar &trans, unsigned faceIdx, unsigned cartElemIdx, const TransMult &transMult) const
Definition: Transmissibility_impl.hpp:1304
Scalar diffusivity(unsigned elemIdx1, unsigned elemIdx2) const
Return the diffusivity for the intersection between two elements.
Definition: Transmissibility_impl.hpp:142
void computeFaceProperties(const Intersection &intersection, const int, const int, const int, const int, DimVector &faceCenterInside, DimVector &faceCenterOutside, DimVector &faceAreaNormal, std::false_type) const
Definition: Transmissibility_impl.hpp:944
void applyEditNncToGridTrans_(const std::unordered_map< std::size_t, int > &globalToLocal)
Multiplies the grid transmissibilities according to EDITNNC.
Definition: Transmissibility_impl.hpp:1100
Dune::FieldMatrix< Scalar, dimWorld, dimWorld > DimMatrix
Definition: Transmissibility.hpp:59
void applyNncMultreg_(const std::unordered_map< std::size_t, int > &globalToLocal)
Definition: Transmissibility_impl.hpp:1202
void applyEditNncrToGridTrans_(const std::unordered_map< std::size_t, int > &globalToLocal)
Resets the grid transmissibilities according to EDITNNCR.
Definition: Transmissibility_impl.hpp:1113
void extractDispersion_()
Definition: Transmissibility_impl.hpp:691
Scalar thermalHalfTrans(unsigned insideElemIdx, unsigned outsideElemIdx) const
Return the thermal "half transmissibility" for the intersection between two elements.
Definition: Transmissibility_impl.hpp:128
void computeHalfDiffusivity_(Scalar &halfDiff, const DimVector &areaNormal, const DimVector &distance, const Scalar &poro) const
Definition: Transmissibility_impl.hpp:1271
void applyAllZMultipliers_(Scalar &trans, unsigned insideFaceIdx, unsigned outsideFaceIdx, unsigned insideCartElemIdx, unsigned outsideCartElemIdx, const TransMult &transMult, const std::array< int, dimWorld > &cartDims, bool pinchTop)
Apply the Multipliers for the case PINCH(4)==TOPBOT.
Definition: Transmissibility_impl.hpp:727
Scalar transmissibilityBoundary(unsigned elemIdx, unsigned boundaryFaceIdx) const
Return the transmissibility for a given boundary segment.
Definition: Transmissibility_impl.hpp:121
void computeHalfTrans_(Scalar &halfTrans, const DimVector &areaNormal, int faceIdx, const DimVector &distance, const DimMatrix &perm) const
Definition: Transmissibility_impl.hpp:1250
std::array< std::vector< double >, 3 > createTransmissibilityArrays_(const std::array< bool, 3 > &is_tran)
Creates TRANS{XYZ} arrays for modification by FieldProps data.
Definition: Transmissibility_impl.hpp:815
void updateFromEclState_(bool global)
Definition: Transmissibility_impl.hpp:778
void update(bool global, const std::function< unsigned int(unsigned int)> &map={}, bool applyNncMultRegT=false)
Compute all transmissibilities.
Definition: Transmissibility_impl.hpp:162
void extractPermeability_()
Definition: Transmissibility_impl.hpp:589
void resetTransmissibilityFromArrays_(const std::array< bool, 3 > &is_tran, const std::array< std::vector< double >, 3 > &trans)
overwrites calculated transmissibilities
Definition: Transmissibility_impl.hpp:882
Dune::FieldVector< Scalar, dimWorld > DimVector
Definition: Transmissibility.hpp:60
Scalar transmissibilityThreshold_
Definition: Transmissibility.hpp:275
void applyEditNncToGridTransHelper_(const std::unordered_map< std::size_t, int > &globalToLocal, const std::string &keyword, const std::vector< NNCdata > &nncs, const std::function< KeywordLocation(const NNCdata &)> &getLocation, const std::function< void(double &, const double &)> &apply)
Definition: Transmissibility_impl.hpp:1126
Scalar transmissibility(unsigned elemIdx1, unsigned elemIdx2) const
Return the transmissibility for the intersection between two elements.
Definition: Transmissibility_impl.hpp:114
const EclipseState & eclState_
Definition: Transmissibility.hpp:270
Scalar dispersivity(unsigned elemIdx1, unsigned elemIdx2) const
Return the dispersivity for the intersection between two elements.
Definition: Transmissibility_impl.hpp:152
Transmissibility(const EclipseState &eclState, const GridView &gridView, const CartesianIndexMapper &cartMapper, const Grid &grid, std::function< std::array< double, dimWorld >(int)> centroids, bool enableEnergy, bool enableDiffusivity, bool enableDispersivity)
Definition: Transmissibility_impl.hpp:89
void applyNncToGridTrans_(const std::unordered_map< std::size_t, int > &cartesianToCompressed)
Definition: Transmissibility_impl.hpp:1032
void applyNtg_(Scalar &trans, unsigned faceIdx, unsigned elemIdx, const std::vector< double > &ntg) const
Definition: Transmissibility_impl.hpp:1338
DimVector distanceVector_(const DimVector ¢er, int faceIdx, unsigned elemIdx, const std::array< std::vector< DimVector >, dimWorld > &axisCentroids) const
Definition: Transmissibility_impl.hpp:1288
void extractPorosity_()
Definition: Transmissibility_impl.hpp:674
void removeNonCartesianTransmissibilities_(bool removeAll)
Definition: Transmissibility_impl.hpp:703
Scalar thermalHalfTransBoundary(unsigned insideElemIdx, unsigned boundaryFaceIdx) const
Definition: Transmissibility_impl.hpp:135
constexpr unsigned elemIdxShift
Definition: Transmissibility_impl.hpp:60
std::pair< std::uint32_t, std::uint32_t > isIdReverse(const std::uint64_t &id)
Definition: Transmissibility_impl.hpp:70
std::uint64_t directionalIsId(std::uint32_t elemIdx1, std::uint32_t elemIdx2)
Definition: Transmissibility_impl.hpp:81
std::uint64_t isId(std::uint32_t elemIdx1, std::uint32_t elemIdx2)
Definition: Transmissibility_impl.hpp:62
Definition: BlackoilPhases.hpp:27