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>
60 constexpr unsigned elemIdxShift = 32;
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);
67 return (elemBIdx<<elemIdxShift) + elemAIdx;
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);
81 std::uint64_t directionalIsId(std::uint32_t elemIdx1, std::uint32_t elemIdx2)
83 return (std::uint64_t(elemIdx1)<<elemIdxShift) + elemIdx2;
87template<
class Gr
id,
class Gr
idView,
class ElementMapper,
class CartesianIndexMapper,
class Scalar>
88Transmissibility<Grid,GridView,ElementMapper,CartesianIndexMapper,Scalar>::
89Transmissibility(
const EclipseState& eclState,
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();
109 transmissibilityThreshold_ = unitSystem.parse(
"Transmissibility").getSIScaling() * 1e-6;
112template<
class Gr
id,
class Gr
idView,
class ElementMapper,
class CartesianIndexMapper,
class Scalar>
116 return trans_.at(details::isId(elemIdx1, elemIdx2));
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>
130 return thermalHalfTrans_.at(details::directionalIsId(insideElemIdx, outsideElemIdx));
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())
147 return diffusivity_.at(details::isId(elemIdx1, elemIdx2));
150template<
class Gr
id,
class Gr
idView,
class ElementMapper,
class CartesianIndexMapper,
class Scalar>
154 if (dispersivity_.empty())
157 return dispersivity_.at(details::isId(elemIdx1, elemIdx2));
160template<
class Gr
id,
class Gr
idView,
class ElementMapper,
class CartesianIndexMapper,
class Scalar>
162update(
bool global,
const TransUpdateQuantities update_quantities,
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;
220 bool pinchOption4ALL;
222 if (comm.rank() == 0) {
223 const auto& eclGrid = eclState_.getInputGrid();
224 pinchActive = eclGrid.isPinchActive();
225 auto pinchTransCalcMode = eclGrid.getPinchOption();
226 useSmallestMultiplier = eclGrid.getMultzOption() == PinchMode::ALL;
227 pinchOption4ALL = (pinchTransCalcMode == PinchMode::ALL);
230 useSmallestMultiplier =
false;
233 if (global && comm.size() > 1) {
234 comm.broadcast(&useSmallestMultiplier, 1, 0);
235 comm.broadcast(&pinchOption4ALL, 1, 0);
236 comm.broadcast(&pinchActive, 1, 0);
240 for (
const auto& elem : elements(gridView_)) {
241 unsigned elemIdx = elemMapper.index(elem);
243 auto isIt = gridView_.ibegin(elem);
244 const auto& isEndIt = gridView_.iend(elem);
245 unsigned boundaryIsIdx = 0;
246 for (; isIt != isEndIt; ++ isIt) {
248 const auto& intersection = *isIt;
251 if (intersection.boundary()) {
253 const auto& geometry = intersection.geometry();
254 const auto& faceCenterInside = geometry.center();
256 auto faceAreaNormal = intersection.centerUnitOuterNormal();
257 faceAreaNormal *= geometry.volume();
259 Scalar transBoundaryIs;
260 computeHalfTrans_(transBoundaryIs,
262 intersection.indexInInside(),
263 distanceVector_(faceCenterInside,
265 permeability_[elemIdx]);
270 unsigned insideCartElemIdx = cartMapper_.cartesianIndex(elemIdx);
271 applyMultipliers_(transBoundaryIs, intersection.indexInInside(), insideCartElemIdx, transMult);
272 transBoundary_[std::make_pair(elemIdx, boundaryIsIdx)] = transBoundaryIs;
276 if (enableEnergy_ && !onlyTrans) {
277 Scalar transBoundaryEnergyIs;
278 computeHalfDiffusivity_(transBoundaryEnergyIs,
280 distanceVector_(faceCenterInside,
283 thermalHalfTransBoundary_[std::make_pair(elemIdx, boundaryIsIdx)] =
284 transBoundaryEnergyIs;
291 if (!intersection.neighbor()) {
298 const auto& outsideElem = intersection.outside();
299 unsigned outsideElemIdx = elemMapper.index(outsideElem);
303 unsigned insideCartElemIdx = this-> lookUpCartesianData_.template getFieldPropCartesianIdx<Grid>(elemIdx);
304 unsigned outsideCartElemIdx = this-> lookUpCartesianData_.template getFieldPropCartesianIdx<Grid>(outsideElemIdx);
312 if (std::tie(insideCartElemIdx, elemIdx) > std::tie(outsideCartElemIdx, outsideElemIdx))
317 int insideFaceIdx = intersection.indexInInside();
318 int outsideFaceIdx = intersection.indexInOutside();
320 if (insideFaceIdx == -1) {
323 assert(outsideFaceIdx == -1);
324 trans_[details::isId(elemIdx, outsideElemIdx)] = 0.0;
325 if (enableEnergy_ && !onlyTrans){
326 thermalHalfTrans_[details::directionalIsId(elemIdx, outsideElemIdx)] = 0.0;
327 thermalHalfTrans_[details::directionalIsId(outsideElemIdx, elemIdx)] = 0.0;
330 if (updateDiffusivity && !onlyTrans) {
331 diffusivity_[details::isId(elemIdx, outsideElemIdx)] = 0.0;
333 if (updateDispersivity && !onlyTrans) {
334 dispersivity_[details::isId(elemIdx, outsideElemIdx)] = 0.0;
339 DimVector faceCenterInside;
340 DimVector faceCenterOutside;
341 DimVector faceAreaNormal;
343 typename std::is_same<Grid, Dune::CpGrid>::type isCpGrid;
344 computeFaceProperties(intersection,
357 computeHalfTrans_(halfTrans1,
360 distanceVector_(faceCenterInside,
362 permeability_[elemIdx]);
363 computeHalfTrans_(halfTrans2,
366 distanceVector_(faceCenterOutside,
368 permeability_[outsideElemIdx]);
370 applyNtg_(halfTrans1, insideFaceIdx, elemIdx, ntg);
371 applyNtg_(halfTrans2, outsideFaceIdx, outsideElemIdx, ntg);
376 if (std::abs(halfTrans1) < 1e-30 || std::abs(halfTrans2) < 1e-30)
380 trans = 1.0 / (1.0/halfTrans1 + 1.0/halfTrans2);
385 if (insideFaceIdx > 3){
386 auto find_layer = [&cartDims](std::size_t cell){
388 auto k = cell / cartDims[1];
391 int kup = find_layer(insideCartElemIdx);
392 int kdown=find_layer(outsideCartElemIdx);
395 assert((kup != kdown) || (insideCartElemIdx == outsideCartElemIdx));
396 if(std::abs(kup -kdown) > 1){
402 if (useSmallestMultiplier)
407 applyAllZMultipliers_(trans, insideFaceIdx, outsideFaceIdx, insideCartElemIdx,
408 outsideCartElemIdx, transMult, cartDims);
412 applyMultipliers_(trans, insideFaceIdx, insideCartElemIdx, transMult);
414 applyMultipliers_(trans, outsideFaceIdx, outsideCartElemIdx, transMult);
418 FaceDir::DirEnum faceDir;
419 switch (insideFaceIdx) {
422 faceDir = FaceDir::XPlus;
427 faceDir = FaceDir::YPlus;
432 faceDir = FaceDir::ZPlus;
436 throw std::logic_error(
"Could not determine a face direction");
439 trans *= transMult.getRegionMultiplier(insideCartElemIdx,
443 trans_[details::isId(elemIdx, outsideElemIdx)] = trans;
446 if (enableEnergy_ && !onlyTrans) {
448 Scalar halfDiffusivity1;
449 Scalar halfDiffusivity2;
451 computeHalfDiffusivity_(halfDiffusivity1,
453 distanceVector_(faceCenterInside,
456 computeHalfDiffusivity_(halfDiffusivity2,
458 distanceVector_(faceCenterOutside,
462 thermalHalfTrans_[details::directionalIsId(elemIdx, outsideElemIdx)] = halfDiffusivity1;
463 thermalHalfTrans_[details::directionalIsId(outsideElemIdx, elemIdx)] = halfDiffusivity2;
467 if (updateDiffusivity && !onlyTrans) {
469 Scalar halfDiffusivity1;
470 Scalar halfDiffusivity2;
472 computeHalfDiffusivity_(halfDiffusivity1,
474 distanceVector_(faceCenterInside,
477 computeHalfDiffusivity_(halfDiffusivity2,
479 distanceVector_(faceCenterOutside,
481 porosity_[outsideElemIdx]);
483 applyNtg_(halfDiffusivity1, insideFaceIdx, elemIdx, ntg);
484 applyNtg_(halfDiffusivity2, outsideFaceIdx, outsideElemIdx, ntg);
488 if (std::abs(halfDiffusivity1) < 1e-30 || std::abs(halfDiffusivity2) < 1e-30)
492 diffusivity = 1.0 / (1.0/halfDiffusivity1 + 1.0/halfDiffusivity2);
495 diffusivity_[details::isId(elemIdx, outsideElemIdx)] = diffusivity;
499 if (updateDispersivity && !onlyTrans) {
501 Scalar halfDispersivity1;
502 Scalar halfDispersivity2;
504 computeHalfDiffusivity_(halfDispersivity1,
506 distanceVector_(faceCenterInside,
508 dispersion_[elemIdx]);
509 computeHalfDiffusivity_(halfDispersivity2,
511 distanceVector_(faceCenterOutside,
513 dispersion_[outsideElemIdx]);
515 applyNtg_(halfDispersivity1, insideFaceIdx, elemIdx, ntg);
516 applyNtg_(halfDispersivity2, outsideFaceIdx, outsideElemIdx, ntg);
520 if (std::abs(halfDispersivity1) < 1e-30 || std::abs(halfDispersivity2) < 1e-30)
524 dispersivity = 1.0 / (1.0/halfDispersivity1 + 1.0/halfDispersivity2);
527 dispersivity_[details::isId(elemIdx, outsideElemIdx)] = dispersivity;
533 this->updateFromEclState_(global);
536 std::unordered_map<std::size_t,int> globalToLocal;
539 for (
const auto& elem : elements(grid_.leafGridView())) {
540 int elemIdx = elemMapper.index(elem);
541 int cartElemIdx = cartMapper_.cartesianIndex(elemIdx);
542 globalToLocal[cartElemIdx] = elemIdx;
551 this->applyEditNncToGridTrans_(globalToLocal);
552 this->applyPinchNncToGridTrans_(globalToLocal);
553 this->applyNncToGridTrans_(globalToLocal);
554 this->applyEditNncrToGridTrans_(globalToLocal);
555 if (applyNncMultregT) {
556 this->applyNncMultreg_(globalToLocal);
558 warnEditNNC_ =
false;
563 this->removeNonCartesianTransmissibilities_(disableNNC);
566template<
class Gr
id,
class Gr
idView,
class ElementMapper,
class CartesianIndexMapper,
class Scalar>
567void Transmissibility<Grid,GridView,ElementMapper,CartesianIndexMapper,Scalar>::
568extractPermeability_()
570 unsigned numElem = gridView_.size(0);
571 permeability_.resize(numElem);
577 const auto& fp = eclState_.fieldProps();
578 if (fp.has_double(
"PERMX")) {
579 const std::vector<double>& permxData = this-> lookUpData_.assignFieldPropsDoubleOnLeaf(fp,
"PERMX");
581 std::vector<double> permyData;
582 if (fp.has_double(
"PERMY"))
583 permyData = this-> lookUpData_.assignFieldPropsDoubleOnLeaf(fp,
"PERMY");
585 permyData = permxData;
587 std::vector<double> permzData;
588 if (fp.has_double(
"PERMZ"))
589 permzData = this-> lookUpData_.assignFieldPropsDoubleOnLeaf(fp,
"PERMZ");
591 permzData = permxData;
593 for (std::size_t dofIdx = 0; dofIdx < numElem; ++ dofIdx) {
594 permeability_[dofIdx] = 0.0;
595 permeability_[dofIdx][0][0] = permxData[dofIdx];
596 permeability_[dofIdx][1][1] = permyData[dofIdx];
597 permeability_[dofIdx][2][2] = permzData[dofIdx];
604 throw std::logic_error(
"Can't read the intrinsic permeability from the ecl state. "
605 "(The PERM{X,Y,Z} keywords are missing)");
608template<
class Gr
id,
class Gr
idView,
class ElementMapper,
class CartesianIndexMapper,
class Scalar>
609void Transmissibility<Grid,GridView,ElementMapper,CartesianIndexMapper,Scalar>::
610extractPermeability_(
const std::function<
unsigned int(
unsigned int)>& map)
612 unsigned numElem = gridView_.size(0);
613 permeability_.resize(numElem);
619 const auto& fp = eclState_.fieldProps();
620 if (fp.has_double(
"PERMX")) {
621 const std::vector<double>& permxData = this-> lookUpData_.assignFieldPropsDoubleOnLeaf(fp,
"PERMX");
623 std::vector<double> permyData;
624 if (fp.has_double(
"PERMY"))
625 permyData = this-> lookUpData_.assignFieldPropsDoubleOnLeaf(fp,
"PERMY");
627 permyData = permxData;
629 std::vector<double> permzData;
630 if (fp.has_double(
"PERMZ"))
631 permzData = this-> lookUpData_.assignFieldPropsDoubleOnLeaf(fp,
"PERMZ");
633 permzData = permxData;
635 for (std::size_t dofIdx = 0; dofIdx < numElem; ++ dofIdx) {
636 permeability_[dofIdx] = 0.0;
637 std::size_t inputDofIdx = map(dofIdx);
638 permeability_[dofIdx][0][0] = permxData[inputDofIdx];
639 permeability_[dofIdx][1][1] = permyData[inputDofIdx];
640 permeability_[dofIdx][2][2] = permzData[inputDofIdx];
647 throw std::logic_error(
"Can't read the intrinsic permeability from the ecl state. "
648 "(The PERM{X,Y,Z} keywords are missing)");
651template<
class Gr
id,
class Gr
idView,
class ElementMapper,
class CartesianIndexMapper,
class Scalar>
652void Transmissibility<Grid,GridView,ElementMapper,CartesianIndexMapper,Scalar>::
659 const auto& fp = eclState_.fieldProps();
660 if (fp.has_double(
"PORO")) {
661 if constexpr (std::is_same_v<Scalar,double>) {
662 porosity_ = this-> lookUpData_.assignFieldPropsDoubleOnLeaf(fp,
"PORO");
664 const auto por = this-> lookUpData_.assignFieldPropsDoubleOnLeaf(fp,
"PORO");
665 porosity_.resize(por.size());
666 std::copy(por.begin(), por.end(), porosity_.begin());
670 throw std::logic_error(
"Can't read the porosityfrom the ecl state. "
671 "(The PORO keywords are missing)");
674template<
class Gr
id,
class Gr
idView,
class ElementMapper,
class CartesianIndexMapper,
class Scalar>
675void Transmissibility<Grid,GridView,ElementMapper,CartesianIndexMapper,Scalar>::
678 if (!enableDispersivity_) {
679 throw std::runtime_error(
"Dispersion disabled at compile time, but the deck "
680 "contains the DISPERC keyword.");
682 const auto& fp = eclState_.fieldProps();
683 if constexpr (std::is_same_v<Scalar,double>) {
684 dispersion_ = this-> lookUpData_.assignFieldPropsDoubleOnLeaf(fp,
"DISPERC");
686 const auto disp = this-> lookUpData_.assignFieldPropsDoubleOnLeaf(fp,
"DISPERC");
687 dispersion_.resize(disp.size());
688 std::copy(disp.begin(), disp.end(), dispersion_.begin());
692template<
class Gr
id,
class Gr
idView,
class ElementMapper,
class CartesianIndexMapper,
class Scalar>
693void Transmissibility<Grid,GridView,ElementMapper,CartesianIndexMapper,Scalar>::
694removeNonCartesianTransmissibilities_(
bool removeAll)
696 const auto& cartDims = cartMapper_.cartesianDimensions();
697 for (
auto&& trans: trans_) {
699 if (removeAll or trans.second < transmissibilityThreshold_) {
700 const auto&
id = trans.first;
701 const auto& elements = details::isIdReverse(
id);
702 int gc1 = std::min(cartMapper_.cartesianIndex(elements.first), cartMapper_.cartesianIndex(elements.second));
703 int gc2 = std::max(cartMapper_.cartesianIndex(elements.first), cartMapper_.cartesianIndex(elements.second));
708 if (gc2 - gc1 == 1 || gc2 - gc1 == cartDims[0] || gc2 - gc1 == cartDims[0]*cartDims[1] || gc2 - gc1 == 0)
716template<
class Gr
id,
class Gr
idView,
class ElementMapper,
class CartesianIndexMapper,
class Scalar>
719 unsigned insideFaceIdx,
720 unsigned outsideFaceIdx,
721 unsigned insideCartElemIdx,
722 unsigned outsideCartElemIdx,
723 const TransMult& transMult,
724 const std::array<int, dimWorld>& cartDims)
726 if(grid_.maxLevel()> 0) {
727 OPM_THROW(std::invalid_argument,
"MULTZ not support with LGRS, yet.");
729 if (insideFaceIdx > 3) {
730 assert(insideFaceIdx==5);
732 assert(outsideCartElemIdx >= insideCartElemIdx);
733 unsigned lastCartElemIdx;
734 if (outsideCartElemIdx == insideCartElemIdx) {
735 lastCartElemIdx = outsideCartElemIdx;
738 lastCartElemIdx = outsideCartElemIdx - cartDims[0]*cartDims[1];
741 Scalar mult = transMult.getMultiplier(lastCartElemIdx , FaceDir::ZPlus) *
742 transMult.getMultiplier(outsideCartElemIdx , FaceDir::ZMinus);
746 for(
auto cartElemIdx = insideCartElemIdx; cartElemIdx < lastCartElemIdx;)
748 auto multiplier = transMult.getMultiplier(cartElemIdx, FaceDir::ZPlus);
749 cartElemIdx += cartDims[0]*cartDims[1];
750 multiplier *= transMult.getMultiplier(cartElemIdx, FaceDir::ZMinus);
751 mult = std::min(mult,
static_cast<Scalar
>(multiplier));
758 applyMultipliers_(trans, insideFaceIdx, insideCartElemIdx, transMult);
759 applyMultipliers_(trans, outsideFaceIdx, outsideCartElemIdx, transMult);
763template<
class Gr
id,
class Gr
idView,
class ElementMapper,
class CartesianIndexMapper,
class Scalar>
767 const FieldPropsManager* fp =
768 (global) ? &(eclState_.fieldProps()) :
769 &(eclState_.globalFieldProps());
771 std::array<bool,3> is_tran {fp->tran_active(
"TRANX"),
772 fp->tran_active(
"TRANY"),
773 fp->tran_active(
"TRANZ")};
775 if( !(is_tran[0] ||is_tran[1] || is_tran[2]) )
781 std::array<std::string, 3> keywords {
"TRANX",
"TRANY",
"TRANZ"};
782 std::array<std::vector<double>,3> trans = createTransmissibilityArrays_(is_tran);
783 auto key = keywords.begin();
784 auto perform = is_tran.begin();
786 for (
auto it = trans.begin(); it != trans.end(); ++it, ++key, ++perform)
789 if(grid_.maxLevel()>0) {
790 OPM_THROW(std::invalid_argument,
"Calculations on TRANX/TRANY/TRANZ arrays are not support with LGRS, yet.");
792 fp->apply_tran(*key, *it);
796 resetTransmissibilityFromArrays_(is_tran, trans);
799template<
class Gr
id,
class Gr
idView,
class ElementMapper,
class CartesianIndexMapper,
class Scalar>
800std::array<std::vector<double>,3>
804 const auto& cartDims = cartMapper_.cartesianDimensions();
805 ElementMapper elemMapper(gridView_, Dune::mcmgElementLayout());
807 auto numElem = gridView_.size(0);
808 std::array<std::vector<double>,3> trans =
809 { std::vector<double>(is_tran[0] ? numElem : 0, 0),
810 std::vector<double>(is_tran[1] ? numElem : 0, 0),
811 std::vector<double>(is_tran[2] ? numElem : 0, 0)};
814 for (
const auto& elem : elements(gridView_)) {
815 for (
const auto& intersection : intersections(gridView_, elem)) {
817 if (!intersection.neighbor())
831 unsigned c1 = elemMapper.index(intersection.inside());
832 unsigned c2 = elemMapper.index(intersection.outside());
833 int gc1 = cartMapper_.cartesianIndex(c1);
834 int gc2 = cartMapper_.cartesianIndex(c2);
835 if (std::tie(gc1, c1) > std::tie(gc2, c2))
842 auto isID = details::isId(c1, c2);
848 if ((gc2 - gc1 == 1 || (gc2 == gc1 && (intersection.indexInInside() == 0 || intersection.indexInInside() == 1)))
849 && cartDims[0] > 1) {
852 trans[0][c1] = trans_[isID];
854 else if ((gc2 - gc1 == cartDims[0] || (gc2 == gc1 && (intersection.indexInInside() == 2 || intersection.indexInInside() == 3)))
855 && cartDims[1] > 1) {
858 trans[1][c1] = trans_[isID];
860 else if (gc2 - gc1 == cartDims[0]*cartDims[1] ||
861 (gc2 == gc1 && (intersection.indexInInside() == 4 || intersection.indexInInside() == 5))) {
864 trans[2][c1] = trans_[isID];
874template<
class Gr
id,
class Gr
idView,
class ElementMapper,
class CartesianIndexMapper,
class Scalar>
877 const std::array<std::vector<double>,3>& trans)
879 const auto& cartDims = cartMapper_.cartesianDimensions();
880 ElementMapper elemMapper(gridView_, Dune::mcmgElementLayout());
883 for (
const auto& elem : elements(gridView_)) {
884 for (
const auto& intersection : intersections(gridView_, elem)) {
885 if (!intersection.neighbor())
899 unsigned c1 = elemMapper.index(intersection.inside());
900 unsigned c2 = elemMapper.index(intersection.outside());
901 int gc1 = cartMapper_.cartesianIndex(c1);
902 int gc2 = cartMapper_.cartesianIndex(c2);
903 if (std::tie(gc1, c1) > std::tie(gc2, c2))
910 auto isID = details::isId(c1, c2);
916 if ((gc2 - gc1 == 1 || (gc2 == gc1 && (intersection.indexInInside() == 0 || intersection.indexInInside() == 1)))
917 && cartDims[0] > 1) {
920 trans_[isID] = trans[0][c1];
922 else if ((gc2 - gc1 == cartDims[0] || (gc2 == gc1 && (intersection.indexInInside() == 2|| intersection.indexInInside() == 3)))
923 && cartDims[1] > 1) {
926 trans_[isID] = trans[1][c1];
928 else if (gc2 - gc1 == cartDims[0]*cartDims[1] ||
929 (gc2 == gc1 && (intersection.indexInInside() == 4 || intersection.indexInInside() == 5))) {
932 trans_[isID] = trans[2][c1];
940template<
class Gr
id,
class Gr
idView,
class ElementMapper,
class CartesianIndexMapper,
class Scalar>
941template<
class Intersection>
948 DimVector& faceCenterInside,
949 DimVector& faceCenterOutside,
950 DimVector& faceAreaNormal,
951 std::false_type)
const
954 const auto& geometry = intersection.geometry();
955 faceCenterInside = geometry.center();
956 faceCenterOutside = faceCenterInside;
958 faceAreaNormal = intersection.centerUnitOuterNormal();
959 faceAreaNormal *= geometry.volume();
962template<
class Gr
id,
class Gr
idView,
class ElementMapper,
class CartesianIndexMapper,
class Scalar>
963template<
class Intersection>
966 const int insideElemIdx,
967 const int insideFaceIdx,
968 const int outsideElemIdx,
969 const int outsideFaceIdx,
970 DimVector& faceCenterInside,
971 DimVector& faceCenterOutside,
972 DimVector& faceAreaNormal,
973 std::true_type)
const
975 int faceIdx = intersection.id();
977 if(grid_.maxLevel() == 0) {
978 faceCenterInside = grid_.faceCenterEcl(insideElemIdx, insideFaceIdx, intersection);
979 faceCenterOutside = grid_.faceCenterEcl(outsideElemIdx, outsideFaceIdx, intersection);
980 faceAreaNormal = grid_.faceAreaNormalEcl(faceIdx);
983 if ((intersection.inside().level() != intersection.outside().level())) {
991 const auto& parentIntersection = grid_.getParentIntersectionFromLgrBoundaryFace(intersection);
992 const auto& parentIntersectionGeometry = parentIntersection.geometry();
996 faceCenterInside = (intersection.inside().level() == 0) ? parentIntersectionGeometry.center() :
997 grid_.faceCenterEcl(insideElemIdx, insideFaceIdx, intersection);
998 faceCenterOutside = (intersection.outside().level() == 0) ? parentIntersectionGeometry.center() :
999 grid_.faceCenterEcl(outsideElemIdx, outsideFaceIdx, intersection);
1007 faceAreaNormal = intersection.centerUnitOuterNormal();
1008 faceAreaNormal *= intersection.geometry().volume();
1011 assert(intersection.inside().level() == intersection.outside().level());
1013 faceCenterInside = grid_.faceCenterEcl(insideElemIdx, insideFaceIdx, intersection);
1014 faceCenterOutside = grid_.faceCenterEcl(outsideElemIdx, outsideFaceIdx, intersection);
1017 if (intersection.inside().level() > 0) {
1018 faceAreaNormal = intersection.centerUnitOuterNormal();
1019 faceAreaNormal *= intersection.geometry().volume();
1022 faceAreaNormal = grid_.faceAreaNormalEcl(faceIdx);
1027template<
class Gr
id,
class Gr
idView,
class ElementMapper,
class CartesianIndexMapper,
class Scalar>
1033 const auto& nnc_input = eclState_.getPinchNNC();
1035 for (
const auto& nncEntry : nnc_input) {
1036 auto c1 = nncEntry.cell1;
1037 auto c2 = nncEntry.cell2;
1038 auto lowIt = cartesianToCompressed.find(c1);
1039 auto highIt = cartesianToCompressed.find(c2);
1040 int low = (lowIt == cartesianToCompressed.end())? -1 : lowIt->second;
1041 int high = (highIt == cartesianToCompressed.end())? -1 : highIt->second;
1044 std::swap(low, high);
1046 if (low == -1 && high == -1)
1050 if (low == -1 || high == -1) {
1058 auto candidate = trans_.find(details::isId(low, high));
1059 if (candidate != trans_.end()) {
1062 candidate->second = nncEntry.trans;
1068template<
class Gr
id,
class Gr
idView,
class ElementMapper,
class CartesianIndexMapper,
class Scalar>
1074 const auto& nnc_input = eclState_.getInputNNC().input();
1076 for (
const auto& nncEntry : nnc_input) {
1077 auto c1 = nncEntry.cell1;
1078 auto c2 = nncEntry.cell2;
1079 auto lowIt = cartesianToCompressed.find(c1);
1080 auto highIt = cartesianToCompressed.find(c2);
1081 int low = (lowIt == cartesianToCompressed.end())? -1 : lowIt->second;
1082 int high = (highIt == cartesianToCompressed.end())? -1 : highIt->second;
1085 std::swap(low, high);
1087 if (low == -1 && high == -1)
1091 if (low == -1 || high == -1) {
1093 std::ostringstream sstr;
1094 sstr <<
"NNC between active and inactive cells ("
1095 << low <<
" -> " << high <<
") with globalcell is (" << c1 <<
"->" << c2 <<
")";
1096 OpmLog::warning(sstr.str());
1101 auto candidate = trans_.find(details::isId(low, high));
1102 if (candidate != trans_.end()) {
1106 candidate->second += nncEntry.trans;
1137template<
class Gr
id,
class Gr
idView,
class ElementMapper,
class CartesianIndexMapper,
class Scalar>
1141 const auto& input = eclState_.getInputNNC();
1142 applyEditNncToGridTransHelper_(globalToLocal,
"EDITNNC",
1144 [&input](
const NNCdata& nnc){
1145 return input.edit_location(nnc);},
1147 [](Scalar& trans,
const Scalar& rhs){ trans *= rhs;});
1150template<
class Gr
id,
class Gr
idView,
class ElementMapper,
class CartesianIndexMapper,
class Scalar>
1154 const auto& input = eclState_.getInputNNC();
1155 applyEditNncToGridTransHelper_(globalToLocal,
"EDITNNCR",
1157 [&input](
const NNCdata& nnc){
1158 return input.editr_location(nnc);},
1160 [](Scalar& trans,
const Scalar& rhs){ trans = rhs;});
1163template<
class Gr
id,
class Gr
idView,
class ElementMapper,
class CartesianIndexMapper,
class Scalar>
1166 const std::string& keyword,
1167 const std::vector<NNCdata>& nncs,
1168 const std::function<KeywordLocation(
const NNCdata&)>& getLocation,
1169 const std::function<
void(Scalar&,
const Scalar&)>& apply)
1173 const auto& cartDims = cartMapper_.cartesianDimensions();
1175 auto format_ijk = [&cartDims](std::size_t cell) -> std::string {
1176 auto i = cell % cartDims[0]; cell /= cartDims[0];
1177 auto j = cell % cartDims[1];
1178 auto k = cell / cartDims[1];
1180 return fmt::format(
"({},{},{})", i + 1,j + 1,k + 1);
1183 auto print_warning = [&format_ijk, &getLocation, &keyword] (
const NNCdata& nnc) {
1184 const auto& location = getLocation( nnc );
1185 auto warning = fmt::format(
"Problem with {} keyword\n"
1187 "No NNC defined for connection {} -> {}", keyword, location.filename,
1188 location.lineno, format_ijk(nnc.cell1), format_ijk(nnc.cell2));
1189 OpmLog::warning(keyword, warning);
1195 auto nnc = nncs.begin();
1196 auto end = nncs.end();
1197 std::size_t warning_count = 0;
1198 while (nnc != end) {
1199 auto c1 = nnc->cell1;
1200 auto c2 = nnc->cell2;
1201 auto lowIt = globalToLocal.find(c1);
1202 auto highIt = globalToLocal.find(c2);
1204 if (lowIt == globalToLocal.end() || highIt == globalToLocal.end()) {
1206 if ( lowIt != highIt && warnEditNNC_) {
1207 print_warning(*nnc);
1214 auto low = lowIt->second, high = highIt->second;
1217 std::swap(low, high);
1219 auto candidate = trans_.find(details::isId(low, high));
1220 if (candidate == trans_.end() && warnEditNNC_) {
1221 print_warning(*nnc);
1227 while (nnc!= end && c1==nnc->cell1 && c2==nnc->cell2) {
1228 apply(candidate->second, nnc->trans);
1234 if (warning_count > 0) {
1235 auto warning = fmt::format(
"Problems with {} keyword\n"
1236 "A total of {} connections not defined in grid", keyword, warning_count);
1237 OpmLog::warning(warning);
1241template<
class Gr
id,
class Gr
idView,
class ElementMapper,
class CartesianIndexMapper,
class Scalar>
1243Transmissibility<Grid,GridView,ElementMapper,CartesianIndexMapper,Scalar>::
1244applyNncMultreg_(
const std::unordered_map<std::size_t,int>& cartesianToCompressed)
1246 const auto& inputNNC = this->eclState_.getInputNNC();
1247 const auto& transMult = this->eclState_.getTransMult();
1249 auto compressedIdx = [&cartesianToCompressed](
const std::size_t globIdx)
1251 auto ixPos = cartesianToCompressed.find(globIdx);
1252 return (ixPos == cartesianToCompressed.end()) ? -1 : ixPos->second;
1266 for (
const auto& nncList : { &NNC::input, &NNC::editr }) {
1267 for (
const auto& nncEntry : (inputNNC.*nncList)()) {
1268 const auto c1 = nncEntry.cell1;
1269 const auto c2 = nncEntry.cell2;
1271 auto low = compressedIdx(c1);
1272 auto high = compressedIdx(c2);
1274 if ((low == -1) || (high == -1)) {
1279 std::swap(low, high);
1282 auto candidate = this->trans_.find(details::isId(low, high));
1283 if (candidate != this->trans_.end()) {
1284 candidate->second *= transMult.getRegionMultiplierNNC(c1, c2);
1290template<
class Gr
id,
class Gr
idView,
class ElementMapper,
class CartesianIndexMapper,
class Scalar>
1291void Transmissibility<Grid,GridView,ElementMapper,CartesianIndexMapper,Scalar>::
1292computeHalfTrans_(Scalar& halfTrans,
1293 const DimVector& areaNormal,
1295 const DimVector& distance,
1296 const DimMatrix& perm)
const
1298 assert(faceIdx >= 0);
1299 unsigned dimIdx = faceIdx/2;
1300 assert(dimIdx < dimWorld);
1301 halfTrans = perm[dimIdx][dimIdx];
1302 halfTrans *= std::abs(Dune::dot(areaNormal, distance));
1303 halfTrans /= distance.two_norm2();
1306template<
class Gr
id,
class Gr
idView,
class ElementMapper,
class CartesianIndexMapper,
class Scalar>
1307void Transmissibility<Grid,GridView,ElementMapper,CartesianIndexMapper,Scalar>::
1308computeHalfDiffusivity_(Scalar& halfDiff,
1309 const DimVector& areaNormal,
1310 const DimVector& distance,
1311 const Scalar& poro)
const
1314 halfDiff *= std::abs(Dune::dot(areaNormal, distance));
1315 halfDiff /= distance.two_norm2();
1318template<
class Gr
id,
class Gr
idView,
class ElementMapper,
class CartesianIndexMapper,
class Scalar>
1319typename Transmissibility<Grid,GridView,ElementMapper,CartesianIndexMapper,Scalar>::DimVector
1320Transmissibility<Grid,GridView,ElementMapper,CartesianIndexMapper,Scalar>::
1321distanceVector_(
const DimVector& faceCenter,
1322 const unsigned& cellIdx)
const
1324 const auto& cellCenter = centroids_(cellIdx);
1325 DimVector x = faceCenter;
1326 for (
unsigned dimIdx = 0; dimIdx < dimWorld; ++dimIdx)
1327 x[dimIdx] -= cellCenter[dimIdx];
1332template<
class Gr
id,
class Gr
idView,
class ElementMapper,
class CartesianIndexMapper,
class Scalar>
1333void Transmissibility<Grid,GridView,ElementMapper,CartesianIndexMapper,Scalar>::
1334applyMultipliers_(Scalar& trans,
1336 unsigned cartElemIdx,
1337 const TransMult& transMult)
const
1344 trans *= transMult.getMultiplier(cartElemIdx, FaceDir::XMinus);
1347 trans *= transMult.getMultiplier(cartElemIdx, FaceDir::XPlus);
1351 trans *= transMult.getMultiplier(cartElemIdx, FaceDir::YMinus);
1354 trans *= transMult.getMultiplier(cartElemIdx, FaceDir::YPlus);
1358 trans *= transMult.getMultiplier(cartElemIdx, FaceDir::ZMinus);
1361 trans *= transMult.getMultiplier(cartElemIdx, FaceDir::ZPlus);
1366template<
class Gr
id,
class Gr
idView,
class ElementMapper,
class CartesianIndexMapper,
class Scalar>
1367void Transmissibility<Grid,GridView,ElementMapper,CartesianIndexMapper,Scalar>::
1368applyNtg_(Scalar& trans,
1371 const std::vector<double>& ntg)
const
1378 trans *= ntg[elemIdx];
1381 trans *= ntg[elemIdx];
1385 trans *= ntg[elemIdx];
1388 trans *= ntg[elemIdx];
Definition Transmissibility.hpp:54
This file contains a set of helper functions used by VFPProd / VFPInj.
Definition blackoilboundaryratevector.hh:37