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>
163 const std::function<
unsigned int(
unsigned int)>& map,
const bool applyNncMultregT)
166 const bool onlyTrans = (update_quantities == TransUpdateQuantities::Trans);
167 const auto& cartDims = cartMapper_.cartesianDimensions();
168 const auto& transMult = eclState_.getTransMult();
169 const auto& comm = gridView_.comm();
170 ElementMapper elemMapper(gridView_, Dune::mcmgElementLayout());
172 unsigned numElements = elemMapper.size();
174 const std::vector<double>& ntg = this->lookUpData_.assignFieldPropsDoubleOnLeaf(eclState_.fieldProps(),
"NTG");
175 const bool updateDiffusivity = eclState_.getSimulationConfig().isDiffusive();
176 const bool updateDispersivity = eclState_.getSimulationConfig().rock_config().dispersion();
178 const bool disableNNC = eclState_.getSimulationConfig().useNONNC();
181 extractPermeability_(map);
183 extractPermeability_();
190 trans_.reserve(numElements*3*1.05);
192 transBoundary_.clear();
195 if (enableEnergy_ && !onlyTrans) {
196 thermalHalfTrans_.clear();
197 thermalHalfTrans_.reserve(numElements*6*1.05);
199 thermalHalfTransBoundary_.clear();
203 if (updateDiffusivity && !onlyTrans) {
204 diffusivity_.clear();
205 diffusivity_.reserve(numElements*3*1.05);
210 if (updateDispersivity && !onlyTrans) {
211 dispersivity_.clear();
212 dispersivity_.reserve(numElements*3*1.05);
213 extractDispersion_();
219 bool useSmallestMultiplier;
221 if (comm.rank() == 0) {
222 const auto& eclGrid = eclState_.getInputGrid();
223 pinchActive = eclGrid.isPinchActive();
224 useSmallestMultiplier = eclGrid.getMultzOption() == PinchMode::ALL;
226 if (global && comm.size() > 1) {
227 comm.broadcast(&useSmallestMultiplier, 1, 0);
228 comm.broadcast(&pinchActive, 1, 0);
232 for (
const auto& elem : elements(gridView_)) {
233 unsigned elemIdx = elemMapper.index(elem);
235 auto isIt = gridView_.ibegin(elem);
236 const auto& isEndIt = gridView_.iend(elem);
237 unsigned boundaryIsIdx = 0;
238 for (; isIt != isEndIt; ++ isIt) {
240 const auto& intersection = *isIt;
243 if (intersection.boundary()) {
245 const auto& geometry = intersection.geometry();
246 const auto& faceCenterInside = geometry.center();
248 auto faceAreaNormal = intersection.centerUnitOuterNormal();
249 faceAreaNormal *= geometry.volume();
251 Scalar transBoundaryIs;
252 computeHalfTrans_(transBoundaryIs,
254 intersection.indexInInside(),
255 distanceVector_(faceCenterInside,
257 permeability_[elemIdx]);
262 unsigned insideCartElemIdx = cartMapper_.cartesianIndex(elemIdx);
263 applyMultipliers_(transBoundaryIs, intersection.indexInInside(), insideCartElemIdx, transMult);
264 transBoundary_[std::make_pair(elemIdx, boundaryIsIdx)] = transBoundaryIs;
268 if (enableEnergy_ && !onlyTrans) {
269 Scalar transBoundaryEnergyIs;
270 computeHalfDiffusivity_(transBoundaryEnergyIs,
272 distanceVector_(faceCenterInside,
275 thermalHalfTransBoundary_[std::make_pair(elemIdx, boundaryIsIdx)] =
276 transBoundaryEnergyIs;
283 if (!intersection.neighbor()) {
290 const auto& outsideElem = intersection.outside();
291 unsigned outsideElemIdx = elemMapper.index(outsideElem);
295 unsigned insideCartElemIdx = this-> lookUpCartesianData_.template getFieldPropCartesianIdx<Grid>(elemIdx);
296 unsigned outsideCartElemIdx = this-> lookUpCartesianData_.template getFieldPropCartesianIdx<Grid>(outsideElemIdx);
304 if (std::tie(insideCartElemIdx, elemIdx) > std::tie(outsideCartElemIdx, outsideElemIdx))
309 int insideFaceIdx = intersection.indexInInside();
310 int outsideFaceIdx = intersection.indexInOutside();
312 if (insideFaceIdx == -1) {
315 assert(outsideFaceIdx == -1);
317 if (enableEnergy_ && !onlyTrans){
322 if (updateDiffusivity && !onlyTrans) {
325 if (updateDispersivity && !onlyTrans) {
335 typename std::is_same<Grid, Dune::CpGrid>::type isCpGrid;
336 computeFaceProperties(intersection,
349 computeHalfTrans_(halfTrans1,
352 distanceVector_(faceCenterInside,
354 permeability_[elemIdx]);
355 computeHalfTrans_(halfTrans2,
358 distanceVector_(faceCenterOutside,
360 permeability_[outsideElemIdx]);
362 applyNtg_(halfTrans1, insideFaceIdx, elemIdx, ntg);
363 applyNtg_(halfTrans2, outsideFaceIdx, outsideElemIdx, ntg);
368 if (std::abs(halfTrans1) < 1e-30 || std::abs(halfTrans2) < 1e-30)
372 trans = 1.0 / (1.0/halfTrans1 + 1.0/halfTrans2);
377 if (insideFaceIdx > 3){
378 auto find_layer = [&cartDims](std::size_t cell){
380 auto k = cell / cartDims[1];
383 int kup = find_layer(insideCartElemIdx);
384 int kdown=find_layer(outsideCartElemIdx);
387 assert((kup != kdown) || (insideCartElemIdx == outsideCartElemIdx));
388 if(std::abs(kup -kdown) > 1){
394 if (useSmallestMultiplier)
398 applyAllZMultipliers_(trans, insideFaceIdx, outsideFaceIdx, insideCartElemIdx,
399 outsideCartElemIdx, transMult, cartDims,
404 applyMultipliers_(trans, insideFaceIdx, insideCartElemIdx, transMult);
406 applyMultipliers_(trans, outsideFaceIdx, outsideCartElemIdx, transMult);
410 FaceDir::DirEnum faceDir;
411 switch (insideFaceIdx) {
414 faceDir = FaceDir::XPlus;
419 faceDir = FaceDir::YPlus;
424 faceDir = FaceDir::ZPlus;
428 throw std::logic_error(
"Could not determine a face direction");
431 trans *= transMult.getRegionMultiplier(insideCartElemIdx,
438 if (enableEnergy_ && !onlyTrans) {
440 Scalar halfDiffusivity1;
441 Scalar halfDiffusivity2;
443 computeHalfDiffusivity_(halfDiffusivity1,
445 distanceVector_(faceCenterInside,
448 computeHalfDiffusivity_(halfDiffusivity2,
450 distanceVector_(faceCenterOutside,
459 if (updateDiffusivity && !onlyTrans) {
461 Scalar halfDiffusivity1;
462 Scalar halfDiffusivity2;
464 computeHalfDiffusivity_(halfDiffusivity1,
466 distanceVector_(faceCenterInside,
469 computeHalfDiffusivity_(halfDiffusivity2,
471 distanceVector_(faceCenterOutside,
473 porosity_[outsideElemIdx]);
475 applyNtg_(halfDiffusivity1, insideFaceIdx, elemIdx, ntg);
476 applyNtg_(halfDiffusivity2, outsideFaceIdx, outsideElemIdx, ntg);
480 if (std::abs(halfDiffusivity1) < 1e-30 || std::abs(halfDiffusivity2) < 1e-30)
484 diffusivity = 1.0 / (1.0/halfDiffusivity1 + 1.0/halfDiffusivity2);
487 diffusivity_[
details::isId(elemIdx, outsideElemIdx)] = diffusivity;
491 if (updateDispersivity && !onlyTrans) {
493 Scalar halfDispersivity1;
494 Scalar halfDispersivity2;
496 computeHalfDiffusivity_(halfDispersivity1,
498 distanceVector_(faceCenterInside,
500 dispersion_[elemIdx]);
501 computeHalfDiffusivity_(halfDispersivity2,
503 distanceVector_(faceCenterOutside,
505 dispersion_[outsideElemIdx]);
507 applyNtg_(halfDispersivity1, insideFaceIdx, elemIdx, ntg);
508 applyNtg_(halfDispersivity2, outsideFaceIdx, outsideElemIdx, ntg);
512 if (std::abs(halfDispersivity1) < 1e-30 || std::abs(halfDispersivity2) < 1e-30)
516 dispersivity = 1.0 / (1.0/halfDispersivity1 + 1.0/halfDispersivity2);
519 dispersivity_[
details::isId(elemIdx, outsideElemIdx)] = dispersivity;
525 this->updateFromEclState_(global);
528 std::unordered_map<std::size_t,int> globalToLocal;
531 for (
const auto& elem : elements(grid_.leafGridView())) {
532 int elemIdx = elemMapper.index(elem);
533 int cartElemIdx = cartMapper_.cartesianIndex(elemIdx);
534 globalToLocal[cartElemIdx] = elemIdx;
543 this->applyEditNncToGridTrans_(globalToLocal);
544 this->applyNncToGridTrans_(globalToLocal);
545 this->applyEditNncrToGridTrans_(globalToLocal);
546 if (applyNncMultregT) {
547 this->applyNncMultreg_(globalToLocal);
549 warnEditNNC_ =
false;
554 this->removeNonCartesianTransmissibilities_(disableNNC);
557template<
class Gr
id,
class Gr
idView,
class ElementMapper,
class CartesianIndexMapper,
class Scalar>
561 unsigned numElem = gridView_.size(0);
562 permeability_.resize(numElem);
568 const auto& fp = eclState_.fieldProps();
569 if (fp.has_double(
"PERMX")) {
570 const std::vector<double>& permxData = this-> lookUpData_.assignFieldPropsDoubleOnLeaf(fp,
"PERMX");
572 std::vector<double> permyData;
573 if (fp.has_double(
"PERMY"))
574 permyData = this-> lookUpData_.assignFieldPropsDoubleOnLeaf(fp,
"PERMY");
576 permyData = permxData;
578 std::vector<double> permzData;
579 if (fp.has_double(
"PERMZ"))
580 permzData = this-> lookUpData_.assignFieldPropsDoubleOnLeaf(fp,
"PERMZ");
582 permzData = permxData;
584 for (std::size_t dofIdx = 0; dofIdx < numElem; ++ dofIdx) {
585 permeability_[dofIdx] = 0.0;
586 permeability_[dofIdx][0][0] = permxData[dofIdx];
587 permeability_[dofIdx][1][1] = permyData[dofIdx];
588 permeability_[dofIdx][2][2] = permzData[dofIdx];
595 throw std::logic_error(
"Can't read the intrinsic permeability from the ecl state. "
596 "(The PERM{X,Y,Z} keywords are missing)");
599template<
class Gr
id,
class Gr
idView,
class ElementMapper,
class CartesianIndexMapper,
class Scalar>
603 unsigned numElem = gridView_.size(0);
604 permeability_.resize(numElem);
610 const auto& fp = eclState_.fieldProps();
611 if (fp.has_double(
"PERMX")) {
612 const std::vector<double>& permxData = this-> lookUpData_.assignFieldPropsDoubleOnLeaf(fp,
"PERMX");
614 std::vector<double> permyData;
615 if (fp.has_double(
"PERMY"))
616 permyData = this-> lookUpData_.assignFieldPropsDoubleOnLeaf(fp,
"PERMY");
618 permyData = permxData;
620 std::vector<double> permzData;
621 if (fp.has_double(
"PERMZ"))
622 permzData = this-> lookUpData_.assignFieldPropsDoubleOnLeaf(fp,
"PERMZ");
624 permzData = permxData;
626 for (std::size_t dofIdx = 0; dofIdx < numElem; ++ dofIdx) {
627 permeability_[dofIdx] = 0.0;
628 std::size_t inputDofIdx = map(dofIdx);
629 permeability_[dofIdx][0][0] = permxData[inputDofIdx];
630 permeability_[dofIdx][1][1] = permyData[inputDofIdx];
631 permeability_[dofIdx][2][2] = permzData[inputDofIdx];
638 throw std::logic_error(
"Can't read the intrinsic permeability from the ecl state. "
639 "(The PERM{X,Y,Z} keywords are missing)");
642template<
class Gr
id,
class Gr
idView,
class ElementMapper,
class CartesianIndexMapper,
class Scalar>
650 const auto& fp = eclState_.fieldProps();
651 if (fp.has_double(
"PORO")) {
652 if constexpr (std::is_same_v<Scalar,double>) {
653 porosity_ = this-> lookUpData_.assignFieldPropsDoubleOnLeaf(fp,
"PORO");
655 const auto por = this-> lookUpData_.assignFieldPropsDoubleOnLeaf(fp,
"PORO");
656 porosity_.resize(por.size());
657 std::copy(por.begin(), por.end(), porosity_.begin());
661 throw std::logic_error(
"Can't read the porosityfrom the ecl state. "
662 "(The PORO keywords are missing)");
665template<
class Gr
id,
class Gr
idView,
class ElementMapper,
class CartesianIndexMapper,
class Scalar>
669 if (!enableDispersivity_) {
670 throw std::runtime_error(
"Dispersion disabled at compile time, but the deck "
671 "contains the DISPERC keyword.");
673 const auto& fp = eclState_.fieldProps();
674 if constexpr (std::is_same_v<Scalar,double>) {
675 dispersion_ = this-> lookUpData_.assignFieldPropsDoubleOnLeaf(fp,
"DISPERC");
677 const auto disp = this-> lookUpData_.assignFieldPropsDoubleOnLeaf(fp,
"DISPERC");
678 dispersion_.resize(disp.size());
679 std::copy(disp.begin(), disp.end(), dispersion_.begin());
683template<
class Gr
id,
class Gr
idView,
class ElementMapper,
class CartesianIndexMapper,
class Scalar>
687 const auto& cartDims = cartMapper_.cartesianDimensions();
688 for (
auto&& trans: trans_) {
690 if (removeAll or trans.second < transmissibilityThreshold_) {
691 const auto&
id = trans.first;
693 int gc1 = std::min(cartMapper_.cartesianIndex(elements.first), cartMapper_.cartesianIndex(elements.second));
694 int gc2 = std::max(cartMapper_.cartesianIndex(elements.first), cartMapper_.cartesianIndex(elements.second));
699 if (gc2 - gc1 == 1 || gc2 - gc1 == cartDims[0] || gc2 - gc1 == cartDims[0]*cartDims[1] || gc2 - gc1 == 0)
707template<
class Gr
id,
class Gr
idView,
class ElementMapper,
class CartesianIndexMapper,
class Scalar>
710 unsigned insideFaceIdx,
711 unsigned outsideFaceIdx,
712 unsigned insideCartElemIdx,
713 unsigned outsideCartElemIdx,
714 const TransMult& transMult,
715 const std::array<int, dimWorld>& cartDims,
718 if(grid_.maxLevel()> 0) {
719 OPM_THROW(std::invalid_argument,
"MULTZ not support with LGRS, yet.");
721 if (insideFaceIdx > 3) {
722 assert(insideFaceIdx==5);
724 assert(outsideCartElemIdx >= insideCartElemIdx);
725 unsigned lastCartElemIdx;
726 if (outsideCartElemIdx == insideCartElemIdx) {
727 lastCartElemIdx = outsideCartElemIdx;
730 lastCartElemIdx = outsideCartElemIdx - cartDims[0]*cartDims[1];
733 Scalar mult = transMult.getMultiplier(lastCartElemIdx , FaceDir::ZPlus) *
734 transMult.getMultiplier(outsideCartElemIdx , FaceDir::ZMinus);
740 for(
auto cartElemIdx = insideCartElemIdx; cartElemIdx < lastCartElemIdx;)
742 auto multiplier = transMult.getMultiplier(cartElemIdx, FaceDir::ZPlus);
743 cartElemIdx += cartDims[0]*cartDims[1];
744 multiplier *= transMult.getMultiplier(cartElemIdx, FaceDir::ZMinus);
745 mult = std::min(mult,
static_cast<Scalar
>(multiplier));
753 applyMultipliers_(trans, insideFaceIdx, insideCartElemIdx, transMult);
754 applyMultipliers_(trans, outsideFaceIdx, outsideCartElemIdx, transMult);
758template<
class Gr
id,
class Gr
idView,
class ElementMapper,
class CartesianIndexMapper,
class Scalar>
762 const FieldPropsManager* fp =
763 (global) ? &(eclState_.fieldProps()) :
764 &(eclState_.globalFieldProps());
766 std::array<bool,3> is_tran {fp->tran_active(
"TRANX"),
767 fp->tran_active(
"TRANY"),
768 fp->tran_active(
"TRANZ")};
770 if( !(is_tran[0] ||is_tran[1] || is_tran[2]) )
776 std::array<std::string, 3> keywords {
"TRANX",
"TRANY",
"TRANZ"};
777 std::array<std::vector<double>,3> trans = createTransmissibilityArrays_(is_tran);
778 auto key = keywords.begin();
779 auto perform = is_tran.begin();
781 for (
auto it = trans.begin(); it != trans.end(); ++it, ++key, ++perform)
784 if(grid_.maxLevel()>0) {
785 OPM_THROW(std::invalid_argument,
"Calculations on TRANX/TRANY/TRANZ arrays are not support with LGRS, yet.");
787 fp->apply_tran(*key, *it);
791 resetTransmissibilityFromArrays_(is_tran, trans);
794template<
class Gr
id,
class Gr
idView,
class ElementMapper,
class CartesianIndexMapper,
class Scalar>
795std::array<std::vector<double>,3>
799 const auto& cartDims = cartMapper_.cartesianDimensions();
800 ElementMapper elemMapper(gridView_, Dune::mcmgElementLayout());
802 auto numElem = gridView_.size(0);
803 std::array<std::vector<double>,3> trans =
804 { std::vector<double>(is_tran[0] ? numElem : 0, 0),
805 std::vector<double>(is_tran[1] ? numElem : 0, 0),
806 std::vector<double>(is_tran[2] ? numElem : 0, 0)};
809 for (
const auto& elem : elements(gridView_)) {
810 for (
const auto& intersection : intersections(gridView_, elem)) {
812 if (!intersection.neighbor())
826 unsigned c1 = elemMapper.index(intersection.inside());
827 unsigned c2 = elemMapper.index(intersection.outside());
828 int gc1 = cartMapper_.cartesianIndex(c1);
829 int gc2 = cartMapper_.cartesianIndex(c2);
830 if (std::tie(gc1, c1) > std::tie(gc2, c2))
843 if ((gc2 - gc1 == 1 || (gc2 == gc1 && (intersection.indexInInside() == 0 || intersection.indexInInside() == 1)))
844 && cartDims[0] > 1) {
847 trans[0][c1] = trans_[isID];
849 else if ((gc2 - gc1 == cartDims[0] || (gc2 == gc1 && (intersection.indexInInside() == 2 || intersection.indexInInside() == 3)))
850 && cartDims[1] > 1) {
853 trans[1][c1] = trans_[isID];
855 else if (gc2 - gc1 == cartDims[0]*cartDims[1] ||
856 (gc2 == gc1 && (intersection.indexInInside() == 4 || intersection.indexInInside() == 5))) {
859 trans[2][c1] = trans_[isID];
869template<
class Gr
id,
class Gr
idView,
class ElementMapper,
class CartesianIndexMapper,
class Scalar>
872 const std::array<std::vector<double>,3>& trans)
874 const auto& cartDims = cartMapper_.cartesianDimensions();
875 ElementMapper elemMapper(gridView_, Dune::mcmgElementLayout());
878 for (
const auto& elem : elements(gridView_)) {
879 for (
const auto& intersection : intersections(gridView_, elem)) {
880 if (!intersection.neighbor())
894 unsigned c1 = elemMapper.index(intersection.inside());
895 unsigned c2 = elemMapper.index(intersection.outside());
896 int gc1 = cartMapper_.cartesianIndex(c1);
897 int gc2 = cartMapper_.cartesianIndex(c2);
898 if (std::tie(gc1, c1) > std::tie(gc2, c2))
911 if ((gc2 - gc1 == 1 || (gc2 == gc1 && (intersection.indexInInside() == 0 || intersection.indexInInside() == 1)))
912 && cartDims[0] > 1) {
915 trans_[isID] = trans[0][c1];
917 else if ((gc2 - gc1 == cartDims[0] || (gc2 == gc1 && (intersection.indexInInside() == 2|| intersection.indexInInside() == 3)))
918 && cartDims[1] > 1) {
921 trans_[isID] = trans[1][c1];
923 else if (gc2 - gc1 == cartDims[0]*cartDims[1] ||
924 (gc2 == gc1 && (intersection.indexInInside() == 4 || intersection.indexInInside() == 5))) {
927 trans_[isID] = trans[2][c1];
935template<
class Gr
id,
class Gr
idView,
class ElementMapper,
class CartesianIndexMapper,
class Scalar>
936template<
class Intersection>
946 std::false_type)
const
949 const auto& geometry = intersection.geometry();
950 faceCenterInside = geometry.center();
951 faceCenterOutside = faceCenterInside;
953 faceAreaNormal = intersection.centerUnitOuterNormal();
954 faceAreaNormal *= geometry.volume();
957template<
class Gr
id,
class Gr
idView,
class ElementMapper,
class CartesianIndexMapper,
class Scalar>
958template<
class Intersection>
961 const int insideElemIdx,
962 const int insideFaceIdx,
963 const int outsideElemIdx,
964 const int outsideFaceIdx,
968 std::true_type)
const
970 int faceIdx = intersection.id();
972 if(grid_.maxLevel() == 0) {
973 faceCenterInside = grid_.faceCenterEcl(insideElemIdx, insideFaceIdx, intersection);
974 faceCenterOutside = grid_.faceCenterEcl(outsideElemIdx, outsideFaceIdx, intersection);
975 faceAreaNormal = grid_.faceAreaNormalEcl(faceIdx);
978 if ((intersection.inside().level() != intersection.outside().level())) {
986 const auto& parentIntersection = grid_.getParentIntersectionFromLgrBoundaryFace(intersection);
987 const auto& parentIntersectionGeometry = parentIntersection.geometry();
991 faceCenterInside = (intersection.inside().level() == 0) ? parentIntersectionGeometry.center() :
992 grid_.faceCenterEcl(insideElemIdx, insideFaceIdx, intersection);
993 faceCenterOutside = (intersection.outside().level() == 0) ? parentIntersectionGeometry.center() :
994 grid_.faceCenterEcl(outsideElemIdx, outsideFaceIdx, intersection);
1002 faceAreaNormal = intersection.centerUnitOuterNormal();
1003 faceAreaNormal *= intersection.geometry().volume();
1006 assert(intersection.inside().level() == intersection.outside().level());
1008 faceCenterInside = grid_.faceCenterEcl(insideElemIdx, insideFaceIdx, intersection);
1009 faceCenterOutside = grid_.faceCenterEcl(outsideElemIdx, outsideFaceIdx, intersection);
1012 if (intersection.inside().level() > 0) {
1013 faceAreaNormal = intersection.centerUnitOuterNormal();
1014 faceAreaNormal *= intersection.geometry().volume();
1017 faceAreaNormal = grid_.faceAreaNormalEcl(faceIdx);
1023template<
class Gr
id,
class Gr
idView,
class ElementMapper,
class CartesianIndexMapper,
class Scalar>
1029 const auto& nnc_input = eclState_.getInputNNC().input();
1031 for (
const auto& nncEntry : nnc_input) {
1032 auto c1 = nncEntry.cell1;
1033 auto c2 = nncEntry.cell2;
1034 auto lowIt = cartesianToCompressed.find(c1);
1035 auto highIt = cartesianToCompressed.find(c2);
1036 int low = (lowIt == cartesianToCompressed.end())? -1 : lowIt->second;
1037 int high = (highIt == cartesianToCompressed.end())? -1 : highIt->second;
1040 std::swap(low, high);
1042 if (low == -1 && high == -1)
1046 if (low == -1 || high == -1) {
1048 std::ostringstream sstr;
1049 sstr <<
"NNC between active and inactive cells ("
1050 << low <<
" -> " << high <<
") with globalcell is (" << c1 <<
"->" << c2 <<
")";
1051 OpmLog::warning(sstr.str());
1057 if (candidate != trans_.end()) {
1061 candidate->second += nncEntry.trans;
1092template<
class Gr
id,
class Gr
idView,
class ElementMapper,
class CartesianIndexMapper,
class Scalar>
1096 const auto& input = eclState_.getInputNNC();
1097 applyEditNncToGridTransHelper_(globalToLocal,
"EDITNNC",
1099 [&input](
const NNCdata& nnc){
1100 return input.edit_location(nnc);},
1102 [](Scalar& trans,
const Scalar& rhs){ trans *= rhs;});
1105template<
class Gr
id,
class Gr
idView,
class ElementMapper,
class CartesianIndexMapper,
class Scalar>
1109 const auto& input = eclState_.getInputNNC();
1110 applyEditNncToGridTransHelper_(globalToLocal,
"EDITNNCR",
1112 [&input](
const NNCdata& nnc){
1113 return input.editr_location(nnc);},
1115 [](Scalar& trans,
const Scalar& rhs){ trans = rhs;});
1118template<
class Gr
id,
class Gr
idView,
class ElementMapper,
class CartesianIndexMapper,
class Scalar>
1121 const std::string& keyword,
1122 const std::vector<NNCdata>& nncs,
1123 const std::function<KeywordLocation(
const NNCdata&)>& getLocation,
1124 const std::function<
void(Scalar&,
const Scalar&)>& apply)
1128 const auto& cartDims = cartMapper_.cartesianDimensions();
1130 auto format_ijk = [&cartDims](std::size_t cell) -> std::string {
1131 auto i = cell % cartDims[0]; cell /= cartDims[0];
1132 auto j = cell % cartDims[1];
1133 auto k = cell / cartDims[1];
1135 return fmt::format(
"({},{},{})", i + 1,j + 1,k + 1);
1138 auto print_warning = [&format_ijk, &getLocation, &keyword] (
const NNCdata& nnc) {
1139 const auto& location = getLocation( nnc );
1140 auto warning = fmt::format(
"Problem with {} keyword\n"
1142 "No NNC defined for connection {} -> {}", keyword, location.filename,
1143 location.lineno, format_ijk(nnc.cell1), format_ijk(nnc.cell2));
1144 OpmLog::warning(keyword, warning);
1150 auto nnc = nncs.begin();
1151 auto end = nncs.end();
1152 std::size_t warning_count = 0;
1153 while (nnc != end) {
1154 auto c1 = nnc->cell1;
1155 auto c2 = nnc->cell2;
1156 auto lowIt = globalToLocal.find(c1);
1157 auto highIt = globalToLocal.find(c2);
1159 if (lowIt == globalToLocal.end() || highIt == globalToLocal.end()) {
1161 if ( lowIt != highIt && warnEditNNC_) {
1162 print_warning(*nnc);
1169 auto low = lowIt->second, high = highIt->second;
1172 std::swap(low, high);
1175 if (candidate == trans_.end() && warnEditNNC_) {
1176 print_warning(*nnc);
1182 while (nnc!= end && c1==nnc->cell1 && c2==nnc->cell2) {
1183 apply(candidate->second, nnc->trans);
1189 if (warning_count > 0) {
1190 auto warning = fmt::format(
"Problems with {} keyword\n"
1191 "A total of {} connections not defined in grid", keyword, warning_count);
1192 OpmLog::warning(warning);
1196template<
class Gr
id,
class Gr
idView,
class ElementMapper,
class CartesianIndexMapper,
class Scalar>
1199applyNncMultreg_(
const std::unordered_map<std::size_t,int>& cartesianToCompressed)
1201 const auto& inputNNC = this->eclState_.getInputNNC();
1202 const auto& transMult = this->eclState_.getTransMult();
1204 auto compressedIdx = [&cartesianToCompressed](
const std::size_t globIdx)
1206 auto ixPos = cartesianToCompressed.find(globIdx);
1207 return (ixPos == cartesianToCompressed.end()) ? -1 : ixPos->second;
1221 for (
const auto& nncList : { &NNC::input, &NNC::editr }) {
1222 for (
const auto& nncEntry : (inputNNC.*nncList)()) {
1223 const auto c1 = nncEntry.cell1;
1224 const auto c2 = nncEntry.cell2;
1226 auto low = compressedIdx(c1);
1227 auto high = compressedIdx(c2);
1229 if ((low == -1) || (high == -1)) {
1234 std::swap(low, high);
1237 auto candidate = this->trans_.find(
details::isId(low, high));
1238 if (candidate != this->trans_.end()) {
1239 candidate->second *= transMult.getRegionMultiplierNNC(c1, c2);
1245template<
class Gr
id,
class Gr
idView,
class ElementMapper,
class CartesianIndexMapper,
class Scalar>
1253 assert(faceIdx >= 0);
1254 unsigned dimIdx = faceIdx/2;
1255 assert(dimIdx < dimWorld);
1256 halfTrans = perm[dimIdx][dimIdx];
1257 halfTrans *= std::abs(Dune::dot(areaNormal, distance));
1258 halfTrans /= distance.two_norm2();
1261template<
class Gr
id,
class Gr
idView,
class ElementMapper,
class CartesianIndexMapper,
class Scalar>
1266 const Scalar& poro)
const
1269 halfDiff *= std::abs(Dune::dot(areaNormal, distance));
1270 halfDiff /= distance.two_norm2();
1273template<
class Gr
id,
class Gr
idView,
class ElementMapper,
class CartesianIndexMapper,
class Scalar>
1277 const unsigned& cellIdx)
const
1279 const auto& cellCenter = centroids_(cellIdx);
1281 for (
unsigned dimIdx = 0; dimIdx < dimWorld; ++dimIdx)
1282 x[dimIdx] -= cellCenter[dimIdx];
1287template<
class Gr
id,
class Gr
idView,
class ElementMapper,
class CartesianIndexMapper,
class Scalar>
1291 unsigned cartElemIdx,
1292 const TransMult& transMult)
const
1299 trans *= transMult.getMultiplier(cartElemIdx, FaceDir::XMinus);
1302 trans *= transMult.getMultiplier(cartElemIdx, FaceDir::XPlus);
1306 trans *= transMult.getMultiplier(cartElemIdx, FaceDir::YMinus);
1309 trans *= transMult.getMultiplier(cartElemIdx, FaceDir::YPlus);
1313 trans *= transMult.getMultiplier(cartElemIdx, FaceDir::ZMinus);
1316 trans *= transMult.getMultiplier(cartElemIdx, FaceDir::ZPlus);
1321template<
class Gr
id,
class Gr
idView,
class ElementMapper,
class CartesianIndexMapper,
class Scalar>
1326 const std::vector<double>& ntg)
const
1333 trans *= ntg[elemIdx];
1336 trans *= ntg[elemIdx];
1340 trans *= ntg[elemIdx];
1343 trans *= ntg[elemIdx];
void applyMultipliers_(Scalar &trans, unsigned faceIdx, unsigned cartElemIdx, const TransMult &transMult) const
Definition: Transmissibility_impl.hpp:1289
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:938
void applyEditNncToGridTrans_(const std::unordered_map< std::size_t, int > &globalToLocal)
Multiplies the grid transmissibilities according to EDITNNC.
Definition: Transmissibility_impl.hpp:1094
Dune::FieldMatrix< Scalar, dimWorld, dimWorld > DimMatrix
Definition: Transmissibility.hpp:59
void update(bool global, TransUpdateQuantities update_quantities=TransUpdateQuantities::All, const std::function< unsigned int(unsigned int)> &map={}, bool applyNncMultRegT=false)
Definition: Transmissibility_impl.hpp:162
DimVector distanceVector_(const DimVector &faceCenter, const unsigned &cellIdx) const
Definition: Transmissibility_impl.hpp:1276
void applyNncMultreg_(const std::unordered_map< std::size_t, int > &globalToLocal)
Definition: Transmissibility_impl.hpp:1199
void applyEditNncrToGridTrans_(const std::unordered_map< std::size_t, int > &globalToLocal)
Resets the grid transmissibilities according to EDITNNCR.
Definition: Transmissibility_impl.hpp:1107
void extractDispersion_()
Definition: Transmissibility_impl.hpp:667
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:1263
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:709
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:1247
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:797
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(Scalar &, const Scalar &)> &apply)
Definition: Transmissibility_impl.hpp:1120
void updateFromEclState_(bool global)
Definition: Transmissibility_impl.hpp:760
void extractPermeability_()
Definition: Transmissibility_impl.hpp:559
TransUpdateQuantities
Compute all transmissibilities.
Definition: Transmissibility.hpp:157
void resetTransmissibilityFromArrays_(const std::array< bool, 3 > &is_tran, const std::array< std::vector< double >, 3 > &trans)
overwrites calculated transmissibilities
Definition: Transmissibility_impl.hpp:871
Dune::FieldVector< Scalar, dimWorld > DimVector
Definition: Transmissibility.hpp:60
Scalar transmissibilityThreshold_
Definition: Transmissibility.hpp:282
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:277
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:1026
void applyNtg_(Scalar &trans, unsigned faceIdx, unsigned elemIdx, const std::vector< double > &ntg) const
Definition: Transmissibility_impl.hpp:1323
void extractPorosity_()
Definition: Transmissibility_impl.hpp:644
void removeNonCartesianTransmissibilities_(bool removeAll)
Definition: Transmissibility_impl.hpp:685
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: blackoilboundaryratevector.hh:37