My Project
Loading...
Searching...
No Matches
Transmissibility_impl.hpp
1// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
2// vi: set et ts=4 sw=4 sts=4:
3/*
4 This file is part of the Open Porous Media project (OPM).
5
6 OPM is free software: you can redistribute it and/or modify
7 it under the terms of the GNU General Public License as published by
8 the Free Software Foundation, either version 2 of the License, or
9 (at your option) any later version.
10
11 OPM is distributed in the hope that it will be useful,
12 but WITHOUT ANY WARRANTY; without even the implied warranty of
13 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 GNU General Public License for more details.
15
16 You should have received a copy of the GNU General Public License
17 along with OPM. If not, see <http://www.gnu.org/licenses/>.
18
19 Consult the COPYING file in the top-level source directory of this
20 module for the precise wording of the license and the list of
21 copyright holders.
22*/
23#ifndef OPM_TRANSMISSIBILITY_IMPL_HPP
24#define OPM_TRANSMISSIBILITY_IMPL_HPP
25
26#include <dune/common/version.hh>
27#include <dune/grid/common/mcmgmapper.hh>
28
29#include <opm/grid/CpGrid.hpp>
30
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>
37
39
40#include <fmt/format.h>
41
42#include <algorithm>
43#include <array>
44#include <cassert>
45#include <cmath>
46#include <cstddef>
47#include <cstdint>
48#include <functional>
49#include <initializer_list>
50#include <sstream>
51#include <stdexcept>
52#include <type_traits>
53#include <utility>
54#include <vector>
55
56namespace Opm {
57
58namespace details {
59
60 constexpr unsigned elemIdxShift = 32; // bits
61
62 std::uint64_t isId(std::uint32_t elemIdx1, std::uint32_t elemIdx2)
63 {
64 std::uint32_t elemAIdx = std::min(elemIdx1, elemIdx2);
65 std::uint64_t elemBIdx = std::max(elemIdx1, elemIdx2);
66
67 return (elemBIdx<<elemIdxShift) + elemAIdx;
68 }
69
70 std::pair<std::uint32_t, std::uint32_t> isIdReverse(const std::uint64_t& id)
71 {
72 // Assigning an unsigned integer to a narrower type discards the most significant bits.
73 // See "The C programming language", section A.6.2.
74 // NOTE that the ordering of element A and B may have changed
75 std::uint32_t elemAIdx = id;
76 std::uint32_t elemBIdx = (id - elemAIdx) >> elemIdxShift;
77
78 return std::make_pair(elemAIdx, elemBIdx);
79 }
80
81 std::uint64_t directionalIsId(std::uint32_t elemIdx1, std::uint32_t elemIdx2)
82 {
83 return (std::uint64_t(elemIdx1)<<elemIdxShift) + elemIdx2;
84 }
85}
86
87template<class Grid, class GridView, class ElementMapper, class CartesianIndexMapper, class Scalar>
88Transmissibility<Grid,GridView,ElementMapper,CartesianIndexMapper,Scalar>::
89Transmissibility(const EclipseState& eclState,
90 const GridView& gridView,
91 const CartesianIndexMapper& cartMapper,
92 const Grid& grid,
93 std::function<std::array<double,dimWorld>(int)> centroids,
94 bool enableEnergy,
97 : eclState_(eclState)
98 , gridView_(gridView)
99 , cartMapper_(cartMapper)
100 , grid_(grid)
101 , centroids_(centroids)
102 , enableEnergy_(enableEnergy)
103 , enableDiffusivity_(enableDiffusivity)
104 , enableDispersivity_(enableDispersivity)
105 , lookUpData_(gridView)
106 , lookUpCartesianData_(gridView, cartMapper)
107{
108 const UnitSystem& unitSystem = eclState_.getDeckUnitSystem();
109 transmissibilityThreshold_ = unitSystem.parse("Transmissibility").getSIScaling() * 1e-6;
110}
111
112template<class Grid, class GridView, class ElementMapper, class CartesianIndexMapper, class Scalar>
114transmissibility(unsigned elemIdx1, unsigned elemIdx2) const
115{
116 return trans_.at(details::isId(elemIdx1, elemIdx2));
117}
118
119template<class Grid, class GridView, class ElementMapper, class CartesianIndexMapper, class Scalar>
121transmissibilityBoundary(unsigned elemIdx, unsigned boundaryFaceIdx) const
122{
123 return transBoundary_.at(std::make_pair(elemIdx, boundaryFaceIdx));
124}
125
126template<class Grid, class GridView, class ElementMapper, class CartesianIndexMapper, class Scalar>
128thermalHalfTrans(unsigned insideElemIdx, unsigned outsideElemIdx) const
129{
130 return thermalHalfTrans_.at(details::directionalIsId(insideElemIdx, outsideElemIdx));
131}
132
133template<class Grid, class GridView, class ElementMapper, class CartesianIndexMapper, class Scalar>
135thermalHalfTransBoundary(unsigned insideElemIdx, unsigned boundaryFaceIdx) const
136{
137 return thermalHalfTransBoundary_.at(std::make_pair(insideElemIdx, boundaryFaceIdx));
138}
139
140template<class Grid, class GridView, class ElementMapper, class CartesianIndexMapper, class Scalar>
142diffusivity(unsigned elemIdx1, unsigned elemIdx2) const
143{
144 if (diffusivity_.empty())
145 return 0.0;
146
147 return diffusivity_.at(details::isId(elemIdx1, elemIdx2));
148}
149
150template<class Grid, class GridView, class ElementMapper, class CartesianIndexMapper, class Scalar>
152dispersivity(unsigned elemIdx1, unsigned elemIdx2) const
153{
154 if (dispersivity_.empty())
155 return 0.0;
156
157 return dispersivity_.at(details::isId(elemIdx1, elemIdx2));
158}
159
160template<class Grid, class GridView, 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)
164{
165 // whether only update the permeability related transmissibility
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());
171
172 unsigned numElements = elemMapper.size();
173 // get the ntg values, the ntg values are modified for the cells merged with minpv
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();
177
178 const bool disableNNC = eclState_.getSimulationConfig().useNONNC();
179
180 if (map)
181 extractPermeability_(map);
182 else
183 extractPermeability_();
184
185 // reserving some space in the hashmap upfront saves quite a bit of time because
186 // resizes are costly for hashmaps and there would be quite a few of them if we
187 // would not have a rough idea of how large the final map will be (the rough idea
188 // is a conforming Cartesian grid).
189 trans_.clear();
190 trans_.reserve(numElements*3*1.05);
191
192 transBoundary_.clear();
193
194 // if energy is enabled, let's do the same for the "thermal half transmissibilities"
195 if (enableEnergy_ && !onlyTrans) {
196 thermalHalfTrans_.clear();
197 thermalHalfTrans_.reserve(numElements*6*1.05);
198
199 thermalHalfTransBoundary_.clear();
200 }
201
202 // if diffusion is enabled, let's do the same for the "diffusivity"
203 if (updateDiffusivity && !onlyTrans) {
204 diffusivity_.clear();
205 diffusivity_.reserve(numElements*3*1.05);
206 extractPorosity_();
207 }
208
209 // if dispersion is enabled, let's do the same for the "dispersivity"
210 if (updateDispersivity && !onlyTrans) {
211 dispersivity_.clear();
212 dispersivity_.reserve(numElements*3*1.05);
213 extractDispersion_();
214 }
215
216 // The MULTZ needs special case if the option is ALL
217 // Then the smallest multiplier is applied.
218 // Default is to apply the top and bottom multiplier
219 bool useSmallestMultiplier;
220 bool pinchOption4ALL;
221 bool pinchActive;
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);
228 if (pinchOption4ALL)
229 {
230 useSmallestMultiplier = false;
231 }
232 }
233 if (global && comm.size() > 1) {
234 comm.broadcast(&useSmallestMultiplier, 1, 0);
235 comm.broadcast(&pinchOption4ALL, 1, 0);
236 comm.broadcast(&pinchActive, 1, 0);
237 }
238
239 // compute the transmissibilities for all intersections
240 for (const auto& elem : elements(gridView_)) {
241 unsigned elemIdx = elemMapper.index(elem);
242
243 auto isIt = gridView_.ibegin(elem);
244 const auto& isEndIt = gridView_.iend(elem);
245 unsigned boundaryIsIdx = 0;
246 for (; isIt != isEndIt; ++ isIt) {
247 // store intersection, this might be costly
248 const auto& intersection = *isIt;
249
250 // deal with grid boundaries
251 if (intersection.boundary()) {
252 // compute the transmissibilty for the boundary intersection
253 const auto& geometry = intersection.geometry();
254 const auto& faceCenterInside = geometry.center();
255
256 auto faceAreaNormal = intersection.centerUnitOuterNormal();
257 faceAreaNormal *= geometry.volume();
258
259 Scalar transBoundaryIs;
260 computeHalfTrans_(transBoundaryIs,
261 faceAreaNormal,
262 intersection.indexInInside(),
263 distanceVector_(faceCenterInside,
264 elemIdx),
265 permeability_[elemIdx]);
266
267 // normally there would be two half-transmissibilities that would be
268 // averaged. on the grid boundary there only is the half
269 // transmissibility of the interior element.
270 unsigned insideCartElemIdx = cartMapper_.cartesianIndex(elemIdx);
271 applyMultipliers_(transBoundaryIs, intersection.indexInInside(), insideCartElemIdx, transMult);
272 transBoundary_[std::make_pair(elemIdx, boundaryIsIdx)] = transBoundaryIs;
273
274 // for boundary intersections we also need to compute the thermal
275 // half transmissibilities
276 if (enableEnergy_ && !onlyTrans) {
277 Scalar transBoundaryEnergyIs;
278 computeHalfDiffusivity_(transBoundaryEnergyIs,
279 faceAreaNormal,
280 distanceVector_(faceCenterInside,
281 elemIdx),
282 1.0);
283 thermalHalfTransBoundary_[std::make_pair(elemIdx, boundaryIsIdx)] =
284 transBoundaryEnergyIs;
285 }
286
287 ++ boundaryIsIdx;
288 continue;
289 }
290
291 if (!intersection.neighbor()) {
292 // elements can be on process boundaries, i.e. they are not on the
293 // domain boundary yet they don't have neighbors.
294 ++ boundaryIsIdx;
295 continue;
296 }
297
298 const auto& outsideElem = intersection.outside();
299 unsigned outsideElemIdx = elemMapper.index(outsideElem);
300
301 // Get the Cartesian indices of the origen cells (parent or equivalent cell on level zero), for CpGrid with LGRs.
302 // For genral grids and no LGRs, get the usual Cartesian Index.
303 unsigned insideCartElemIdx = this-> lookUpCartesianData_.template getFieldPropCartesianIdx<Grid>(elemIdx);
304 unsigned outsideCartElemIdx = this-> lookUpCartesianData_.template getFieldPropCartesianIdx<Grid>(outsideElemIdx);
305
306 // we only need to calculate a face's transmissibility
307 // once...
308 // In a parallel run insideCartElemIdx>outsideCartElemIdx does not imply elemIdx>outsideElemIdx for
309 // ghost cells and we need to use the cartesian index as this will be used when applying Z multipliers
310 // To cover the case where both cells are part of an LGR and as a consequence might have
311 // the same cartesian index, we tie their Cartesian indices and the ones on the leaf grid view.
312 if (std::tie(insideCartElemIdx, elemIdx) > std::tie(outsideCartElemIdx, outsideElemIdx))
313 continue;
314
315 // local indices of the faces of the inside and
316 // outside elements which contain the intersection
317 int insideFaceIdx = intersection.indexInInside();
318 int outsideFaceIdx = intersection.indexInOutside();
319
320 if (insideFaceIdx == -1) {
321 // NNC. Set zero transmissibility, as it will be
322 // *added to* by applyNncToGridTrans_() later.
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;
328 }
329
330 if (updateDiffusivity && !onlyTrans) {
331 diffusivity_[details::isId(elemIdx, outsideElemIdx)] = 0.0;
332 }
333 if (updateDispersivity && !onlyTrans) {
334 dispersivity_[details::isId(elemIdx, outsideElemIdx)] = 0.0;
335 }
336 continue;
337 }
338
339 DimVector faceCenterInside;
340 DimVector faceCenterOutside;
341 DimVector faceAreaNormal;
342
343 typename std::is_same<Grid, Dune::CpGrid>::type isCpGrid;
344 computeFaceProperties(intersection,
345 elemIdx,
346 insideFaceIdx,
347 outsideElemIdx,
348 outsideFaceIdx,
349 faceCenterInside,
350 faceCenterOutside,
351 faceAreaNormal,
352 isCpGrid);
353
354 Scalar halfTrans1;
355 Scalar halfTrans2;
356
357 computeHalfTrans_(halfTrans1,
358 faceAreaNormal,
359 insideFaceIdx,
360 distanceVector_(faceCenterInside,
361 elemIdx),
362 permeability_[elemIdx]);
363 computeHalfTrans_(halfTrans2,
364 faceAreaNormal,
365 outsideFaceIdx,
366 distanceVector_(faceCenterOutside,
367 outsideElemIdx),
368 permeability_[outsideElemIdx]);
369
370 applyNtg_(halfTrans1, insideFaceIdx, elemIdx, ntg);
371 applyNtg_(halfTrans2, outsideFaceIdx, outsideElemIdx, ntg);
372
373 // convert half transmissibilities to full face
374 // transmissibilities using the harmonic mean
375 Scalar trans;
376 if (std::abs(halfTrans1) < 1e-30 || std::abs(halfTrans2) < 1e-30)
377 // avoid division by zero
378 trans = 0.0;
379 else
380 trans = 1.0 / (1.0/halfTrans1 + 1.0/halfTrans2);
381
382 // apply the full face transmissibility multipliers
383 // for the inside ...
384 if(!pinchActive){
385 if (insideFaceIdx > 3){// top or bottom
386 auto find_layer = [&cartDims](std::size_t cell){
387 cell /= cartDims[0];
388 auto k = cell / cartDims[1];
389 return k;
390 };
391 int kup = find_layer(insideCartElemIdx);
392 int kdown=find_layer(outsideCartElemIdx);
393 // When a grid is a CpGrid with LGRs, insideCartElemIdx coincides with outsideCartElemIdx
394 // for cells on the leaf with the same parent cell on level zero.
395 assert((kup != kdown) || (insideCartElemIdx == outsideCartElemIdx));
396 if(std::abs(kup -kdown) > 1){
397 trans = 0.0;
398 }
399 }
400 }
401
402 if (useSmallestMultiplier)
403 {
404 // PINCH(4) == TOPBOT is assumed here as we set useSmallestMultipliers
405 // to false if PINCH(4) == ALL holds
406 // In contrast to the name this will also apply
407 applyAllZMultipliers_(trans, insideFaceIdx, outsideFaceIdx, insideCartElemIdx,
408 outsideCartElemIdx, transMult, cartDims);
409 }
410 else
411 {
412 applyMultipliers_(trans, insideFaceIdx, insideCartElemIdx, transMult);
413 // ... and outside elements
414 applyMultipliers_(trans, outsideFaceIdx, outsideCartElemIdx, transMult);
415 }
416
417 // apply the region multipliers (cf. the MULTREGT keyword)
418 FaceDir::DirEnum faceDir;
419 switch (insideFaceIdx) {
420 case 0:
421 case 1:
422 faceDir = FaceDir::XPlus;
423 break;
424
425 case 2:
426 case 3:
427 faceDir = FaceDir::YPlus;
428 break;
429
430 case 4:
431 case 5:
432 faceDir = FaceDir::ZPlus;
433 break;
434
435 default:
436 throw std::logic_error("Could not determine a face direction");
437 }
438
439 trans *= transMult.getRegionMultiplier(insideCartElemIdx,
440 outsideCartElemIdx,
441 faceDir);
442
443 trans_[details::isId(elemIdx, outsideElemIdx)] = trans;
444
445 // update the "thermal half transmissibility" for the intersection
446 if (enableEnergy_ && !onlyTrans) {
447
448 Scalar halfDiffusivity1;
449 Scalar halfDiffusivity2;
450
451 computeHalfDiffusivity_(halfDiffusivity1,
452 faceAreaNormal,
453 distanceVector_(faceCenterInside,
454 elemIdx),
455 1.0);
456 computeHalfDiffusivity_(halfDiffusivity2,
457 faceAreaNormal,
458 distanceVector_(faceCenterOutside,
459 outsideElemIdx),
460 1.0);
461 //TODO Add support for multipliers
462 thermalHalfTrans_[details::directionalIsId(elemIdx, outsideElemIdx)] = halfDiffusivity1;
463 thermalHalfTrans_[details::directionalIsId(outsideElemIdx, elemIdx)] = halfDiffusivity2;
464 }
465
466 // update the "diffusive half transmissibility" for the intersection
467 if (updateDiffusivity && !onlyTrans) {
468
469 Scalar halfDiffusivity1;
470 Scalar halfDiffusivity2;
471
472 computeHalfDiffusivity_(halfDiffusivity1,
473 faceAreaNormal,
474 distanceVector_(faceCenterInside,
475 elemIdx),
476 porosity_[elemIdx]);
477 computeHalfDiffusivity_(halfDiffusivity2,
478 faceAreaNormal,
479 distanceVector_(faceCenterOutside,
480 outsideElemIdx),
481 porosity_[outsideElemIdx]);
482
483 applyNtg_(halfDiffusivity1, insideFaceIdx, elemIdx, ntg);
484 applyNtg_(halfDiffusivity2, outsideFaceIdx, outsideElemIdx, ntg);
485
486 //TODO Add support for multipliers
487 Scalar diffusivity;
488 if (std::abs(halfDiffusivity1) < 1e-30 || std::abs(halfDiffusivity2) < 1e-30)
489 // avoid division by zero
490 diffusivity = 0.0;
491 else
492 diffusivity = 1.0 / (1.0/halfDiffusivity1 + 1.0/halfDiffusivity2);
493
494
495 diffusivity_[details::isId(elemIdx, outsideElemIdx)] = diffusivity;
496 }
497
498 // update the "dispersivity half transmissibility" for the intersection
499 if (updateDispersivity && !onlyTrans) {
500
501 Scalar halfDispersivity1;
502 Scalar halfDispersivity2;
503
504 computeHalfDiffusivity_(halfDispersivity1,
505 faceAreaNormal,
506 distanceVector_(faceCenterInside,
507 elemIdx),
508 dispersion_[elemIdx]);
509 computeHalfDiffusivity_(halfDispersivity2,
510 faceAreaNormal,
511 distanceVector_(faceCenterOutside,
512 outsideElemIdx),
513 dispersion_[outsideElemIdx]);
514
515 applyNtg_(halfDispersivity1, insideFaceIdx, elemIdx, ntg);
516 applyNtg_(halfDispersivity2, outsideFaceIdx, outsideElemIdx, ntg);
517
518 //TODO Add support for multipliers
519 Scalar dispersivity;
520 if (std::abs(halfDispersivity1) < 1e-30 || std::abs(halfDispersivity2) < 1e-30)
521 // avoid division by zero
522 dispersivity = 0.0;
523 else
524 dispersivity = 1.0 / (1.0/halfDispersivity1 + 1.0/halfDispersivity2);
525
526
527 dispersivity_[details::isId(elemIdx, outsideElemIdx)] = dispersivity;
528 }
529 }
530 }
531
532 // Potentially overwrite and/or modify transmissibilities based on input from deck
533 this->updateFromEclState_(global);
534
535 // Create mapping from global to local index
536 std::unordered_map<std::size_t,int> globalToLocal;
537
538 // Loop over all elements (global grid) and store Cartesian index
539 for (const auto& elem : elements(grid_.leafGridView())) {
540 int elemIdx = elemMapper.index(elem);
541 int cartElemIdx = cartMapper_.cartesianIndex(elemIdx);
542 globalToLocal[cartElemIdx] = elemIdx;
543 }
544
545 if (!disableNNC) {
546 // For EDITNNC and EDITNNCR we warn only once
547 // If transmissibility is used for load balancing this will be done
548 // when computing the gobal transmissibilities and all warnings will
549 // be seen in a parallel. Unfortunately, when we do not use transmissibilities
550 // we will only see warnings for the partition of process 0 and also false positives.
551 this->applyEditNncToGridTrans_(globalToLocal);
552 this->applyPinchNncToGridTrans_(globalToLocal);
553 this->applyNncToGridTrans_(globalToLocal);
554 this->applyEditNncrToGridTrans_(globalToLocal);
555 if (applyNncMultregT) {
556 this->applyNncMultreg_(globalToLocal);
557 }
558 warnEditNNC_ = false;
559 }
560
561 // If disableNNC == true, remove all non-neighbouring transmissibilities.
562 // If disableNNC == false, remove very small non-neighbouring transmissibilities.
563 this->removeNonCartesianTransmissibilities_(disableNNC);
564}
565
566template<class Grid, class GridView, class ElementMapper, class CartesianIndexMapper, class Scalar>
567void Transmissibility<Grid,GridView,ElementMapper,CartesianIndexMapper,Scalar>::
568extractPermeability_()
569{
570 unsigned numElem = gridView_.size(/*codim=*/0);
571 permeability_.resize(numElem);
572
573 // read the intrinsic permeabilities from the eclState. Note that all arrays
574 // provided by eclState are one-per-cell of "uncompressed" grid, whereas the
575 // simulation grid might remove a few elements. (e.g. because it is distributed
576 // over several processes.)
577 const auto& fp = eclState_.fieldProps();
578 if (fp.has_double("PERMX")) {
579 const std::vector<double>& permxData = this-> lookUpData_.assignFieldPropsDoubleOnLeaf(fp, "PERMX");
580
581 std::vector<double> permyData;
582 if (fp.has_double("PERMY"))
583 permyData = this-> lookUpData_.assignFieldPropsDoubleOnLeaf(fp,"PERMY");
584 else
585 permyData = permxData;
586
587 std::vector<double> permzData;
588 if (fp.has_double("PERMZ"))
589 permzData = this-> lookUpData_.assignFieldPropsDoubleOnLeaf(fp,"PERMZ");
590 else
591 permzData = permxData;
592
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];
598 }
599
600 // for now we don't care about non-diagonal entries
601
602 }
603 else
604 throw std::logic_error("Can't read the intrinsic permeability from the ecl state. "
605 "(The PERM{X,Y,Z} keywords are missing)");
606}
607
608template<class Grid, class GridView, class ElementMapper, class CartesianIndexMapper, class Scalar>
609void Transmissibility<Grid,GridView,ElementMapper,CartesianIndexMapper,Scalar>::
610extractPermeability_(const std::function<unsigned int(unsigned int)>& map)
611{
612 unsigned numElem = gridView_.size(/*codim=*/0);
613 permeability_.resize(numElem);
614
615 // read the intrinsic permeabilities from the eclState. Note that all arrays
616 // provided by eclState are one-per-cell of "uncompressed" grid, whereas the
617 // simulation grid might remove a few elements. (e.g. because it is distributed
618 // over several processes.)
619 const auto& fp = eclState_.fieldProps();
620 if (fp.has_double("PERMX")) {
621 const std::vector<double>& permxData = this-> lookUpData_.assignFieldPropsDoubleOnLeaf(fp,"PERMX");
622
623 std::vector<double> permyData;
624 if (fp.has_double("PERMY"))
625 permyData = this-> lookUpData_.assignFieldPropsDoubleOnLeaf(fp,"PERMY");
626 else
627 permyData = permxData;
628
629 std::vector<double> permzData;
630 if (fp.has_double("PERMZ"))
631 permzData = this-> lookUpData_.assignFieldPropsDoubleOnLeaf(fp,"PERMZ");
632 else
633 permzData = permxData;
634
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];
641 }
642
643 // for now we don't care about non-diagonal entries
644
645 }
646 else
647 throw std::logic_error("Can't read the intrinsic permeability from the ecl state. "
648 "(The PERM{X,Y,Z} keywords are missing)");
649}
650
651template<class Grid, class GridView, class ElementMapper, class CartesianIndexMapper, class Scalar>
652void Transmissibility<Grid,GridView,ElementMapper,CartesianIndexMapper,Scalar>::
653extractPorosity_()
654{
655 // read the intrinsic porosity from the eclState. Note that all arrays
656 // provided by eclState are one-per-cell of "uncompressed" grid, whereas the
657 // simulation grid might remove a few elements. (e.g. because it is distributed
658 // over several processes.)
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");
663 } else {
664 const auto por = this-> lookUpData_.assignFieldPropsDoubleOnLeaf(fp,"PORO");
665 porosity_.resize(por.size());
666 std::copy(por.begin(), por.end(), porosity_.begin());
667 }
668 }
669 else
670 throw std::logic_error("Can't read the porosityfrom the ecl state. "
671 "(The PORO keywords are missing)");
672}
673
674template<class Grid, class GridView, class ElementMapper, class CartesianIndexMapper, class Scalar>
675void Transmissibility<Grid,GridView,ElementMapper,CartesianIndexMapper,Scalar>::
676extractDispersion_()
677{
678 if (!enableDispersivity_) {
679 throw std::runtime_error("Dispersion disabled at compile time, but the deck "
680 "contains the DISPERC keyword.");
681 }
682 const auto& fp = eclState_.fieldProps();
683 if constexpr (std::is_same_v<Scalar,double>) {
684 dispersion_ = this-> lookUpData_.assignFieldPropsDoubleOnLeaf(fp,"DISPERC");
685 } else {
686 const auto disp = this-> lookUpData_.assignFieldPropsDoubleOnLeaf(fp,"DISPERC");
687 dispersion_.resize(disp.size());
688 std::copy(disp.begin(), disp.end(), dispersion_.begin());
689 }
690}
691
692template<class Grid, class GridView, class ElementMapper, class CartesianIndexMapper, class Scalar>
693void Transmissibility<Grid,GridView,ElementMapper,CartesianIndexMapper,Scalar>::
694removeNonCartesianTransmissibilities_(bool removeAll)
695{
696 const auto& cartDims = cartMapper_.cartesianDimensions();
697 for (auto&& trans: trans_) {
698 //either remove all NNC transmissibilities or those less than the threshold (by default 1e-6 in the deck's unit system)
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));
704
705 // only adjust the NNCs
706 // When LGRs, all neighbors in the LGR are cartesian neighbours on the level grid representing the LGR.
707 // When elements on the leaf grid view have the same parent cell, gc1 and gc2 coincide.
708 if (gc2 - gc1 == 1 || gc2 - gc1 == cartDims[0] || gc2 - gc1 == cartDims[0]*cartDims[1] || gc2 - gc1 == 0)
709 continue;
710
711 trans.second = 0.0;
712 }
713 }
714}
715
716template<class Grid, class GridView, class ElementMapper, class CartesianIndexMapper, class Scalar>
718applyAllZMultipliers_(Scalar& trans,
719 unsigned insideFaceIdx,
720 unsigned outsideFaceIdx,
721 unsigned insideCartElemIdx,
722 unsigned outsideCartElemIdx,
723 const TransMult& transMult,
724 const std::array<int, dimWorld>& cartDims)
725{
726 if(grid_.maxLevel()> 0) {
727 OPM_THROW(std::invalid_argument, "MULTZ not support with LGRS, yet.");
728 }
729 if (insideFaceIdx > 3) { // top or or bottom
730 assert(insideFaceIdx==5); // as insideCartElemIdx < outsideCartElemIdx holds for the Z column
731 // For CpGrid with LGRs, insideCartElemIdx == outsideCartElemIdx when cells on the leaf have the same parent cell on level zero.
732 assert(outsideCartElemIdx >= insideCartElemIdx);
733 unsigned lastCartElemIdx;
734 if (outsideCartElemIdx == insideCartElemIdx) {
735 lastCartElemIdx = outsideCartElemIdx;
736 }
737 else {
738 lastCartElemIdx = outsideCartElemIdx - cartDims[0]*cartDims[1];
739 }
740 // Last multiplier using (Z+)*(Z-)
741 Scalar mult = transMult.getMultiplier(lastCartElemIdx , FaceDir::ZPlus) *
742 transMult.getMultiplier(outsideCartElemIdx , FaceDir::ZMinus);
743
744 // pick the smallest multiplier using (Z+)*(Z-) while looking down
745 // the pillar until reaching the other end of the connection
746 for(auto cartElemIdx = insideCartElemIdx; cartElemIdx < lastCartElemIdx;)
747 {
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));
752 }
753
754 trans *= mult;
755 }
756 else
757 {
758 applyMultipliers_(trans, insideFaceIdx, insideCartElemIdx, transMult);
759 applyMultipliers_(trans, outsideFaceIdx, outsideCartElemIdx, transMult);
760 }
761}
762
763template<class Grid, class GridView, class ElementMapper, class CartesianIndexMapper, class Scalar>
765updateFromEclState_(bool global)
766{
767 const FieldPropsManager* fp =
768 (global) ? &(eclState_.fieldProps()) :
769 &(eclState_.globalFieldProps());
770
771 std::array<bool,3> is_tran {fp->tran_active("TRANX"),
772 fp->tran_active("TRANY"),
773 fp->tran_active("TRANZ")};
774
775 if( !(is_tran[0] ||is_tran[1] || is_tran[2]) )
776 {
777 // Skip unneeded expensive traversals
778 return;
779 }
780
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();
785
786 for (auto it = trans.begin(); it != trans.end(); ++it, ++key, ++perform)
787 {
788 if(*perform) {
789 if(grid_.maxLevel()>0) {
790 OPM_THROW(std::invalid_argument, "Calculations on TRANX/TRANY/TRANZ arrays are not support with LGRS, yet.");
791 }
792 fp->apply_tran(*key, *it);
793 }
794 }
795
796 resetTransmissibilityFromArrays_(is_tran, trans);
797}
798
799template<class Grid, class GridView, class ElementMapper, class CartesianIndexMapper, class Scalar>
800std::array<std::vector<double>,3>
802createTransmissibilityArrays_(const std::array<bool,3>& is_tran)
803{
804 const auto& cartDims = cartMapper_.cartesianDimensions();
805 ElementMapper elemMapper(gridView_, Dune::mcmgElementLayout());
806
807 auto numElem = gridView_.size(/*codim=*/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)};
812
813 // compute the transmissibilities for all intersections
814 for (const auto& elem : elements(gridView_)) {
815 for (const auto& intersection : intersections(gridView_, elem)) {
816 // store intersection, this might be costly
817 if (!intersection.neighbor())
818 continue; // intersection is on the domain boundary
819
820 // In the EclState TRANX[c1] is transmissibility in X+
821 // direction. we only store transmissibilities in the +
822 // direction. Same for Y and Z. Ordering of compressed (c1,c2) and cartesian index
823 // (gc1, gc2) is coherent (c1 < c2 <=> gc1 < gc2) only in a serial run.
824 // In a parallel run this only holds in the interior as elements in the
825 // ghost overlap region might be ordered after the others. Hence we need
826 // to use the cartesian index to select the compressed index where to store
827 // the transmissibility value.
828 // c1 < c2 <=> gc1 < gc2 is no longer true (even in serial) when the grid is a
829 // CpGrid with LGRs. When cells c1 and c2 have the same parent
830 // cell on level zero, then gc1 == gc2.
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))
836 // we only need to handle each connection once, thank you.
837 // We do this when gc1 is smaller than the other to find the
838 // correct place to store in parallel when ghost/overlap elements
839 // are ordered last
840 continue;
841
842 auto isID = details::isId(c1, c2);
843
844 // For CpGrid with LGRs, when leaf grid view cells with indices c1 and c2
845 // have the same parent cell on level zero, then gc2 - gc1 == 0. In that case,
846 // 'intersection.indexInSIde()' needed to be checked to determine the direction, i.e.
847 // add in the if/else-if 'gc2 == gc1 && intersection.indexInInside() == ... '
848 if ((gc2 - gc1 == 1 || (gc2 == gc1 && (intersection.indexInInside() == 0 || intersection.indexInInside() == 1)))
849 && cartDims[0] > 1) {
850 if (is_tran[0])
851 // set simulator internal transmissibilities to values from inputTranx
852 trans[0][c1] = trans_[isID];
853 }
854 else if ((gc2 - gc1 == cartDims[0] || (gc2 == gc1 && (intersection.indexInInside() == 2 || intersection.indexInInside() == 3)))
855 && cartDims[1] > 1) {
856 if (is_tran[1])
857 // set simulator internal transmissibilities to values from inputTrany
858 trans[1][c1] = trans_[isID];
859 }
860 else if (gc2 - gc1 == cartDims[0]*cartDims[1] ||
861 (gc2 == gc1 && (intersection.indexInInside() == 4 || intersection.indexInInside() == 5))) {
862 if (is_tran[2])
863 // set simulator internal transmissibilities to values from inputTranz
864 trans[2][c1] = trans_[isID];
865 }
866
867 //else.. We don't support modification of NNC at the moment.
868 }
869 }
870
871 return trans;
872}
873
874template<class Grid, class GridView, class ElementMapper, class CartesianIndexMapper, class Scalar>
876resetTransmissibilityFromArrays_(const std::array<bool,3>& is_tran,
877 const std::array<std::vector<double>,3>& trans)
878{
879 const auto& cartDims = cartMapper_.cartesianDimensions();
880 ElementMapper elemMapper(gridView_, Dune::mcmgElementLayout());
881
882 // compute the transmissibilities for all intersections
883 for (const auto& elem : elements(gridView_)) {
884 for (const auto& intersection : intersections(gridView_, elem)) {
885 if (!intersection.neighbor())
886 continue; // intersection is on the domain boundary
887
888 // In the EclState TRANX[c1] is transmissibility in X+
889 // direction. we only store transmissibilities in the +
890 // direction. Same for Y and Z. Ordering of compressed (c1,c2) and cartesian index
891 // (gc1, gc2) is coherent (c1 < c2 <=> gc1 < gc2) only in a serial run.
892 // In a parallel run this only holds in the interior as elements in the
893 // ghost overlap region might be ordered after the others. Hence we need
894 // to use the cartesian index to select the compressed index where to store
895 // the transmissibility value.
896 // c1 < c2 <=> gc1 < gc2 is no longer true (even in serial) when the grid is a
897 // CpGrid with LGRs. When cells c1 and c2 have the same parent
898 // cell on level zero, then gc1 == gc2.
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))
904 // we only need to handle each connection once, thank you.
905 // We do this when gc1 is smaller than the other to find the
906 // correct place to read in parallel when ghost/overlap elements
907 // are ordered last
908 continue;
909
910 auto isID = details::isId(c1, c2);
911
912 // For CpGrid with LGRs, when leaf grid view cells with indices c1 and c2
913 // have the same parent cell on level zero, then gc2 - gc1 == 0. In that case,
914 // 'intersection.indexInSIde()' needed to be checked to determine the direction, i.e.
915 // add in the if/else-if 'gc2 == gc1 && intersection.indexInInside() == ... '
916 if ((gc2 - gc1 == 1 || (gc2 == gc1 && (intersection.indexInInside() == 0 || intersection.indexInInside() == 1)))
917 && cartDims[0] > 1) {
918 if (is_tran[0])
919 // set simulator internal transmissibilities to values from inputTranx
920 trans_[isID] = trans[0][c1];
921 }
922 else if ((gc2 - gc1 == cartDims[0] || (gc2 == gc1 && (intersection.indexInInside() == 2|| intersection.indexInInside() == 3)))
923 && cartDims[1] > 1) {
924 if (is_tran[1])
925 // set simulator internal transmissibilities to values from inputTrany
926 trans_[isID] = trans[1][c1];
927 }
928 else if (gc2 - gc1 == cartDims[0]*cartDims[1] ||
929 (gc2 == gc1 && (intersection.indexInInside() == 4 || intersection.indexInInside() == 5))) {
930 if (is_tran[2])
931 // set simulator internal transmissibilities to values from inputTranz
932 trans_[isID] = trans[2][c1];
933 }
934
935 //else.. We don't support modification of NNC at the moment.
936 }
937 }
938}
939
940template<class Grid, class GridView, class ElementMapper, class CartesianIndexMapper, class Scalar>
941template<class Intersection>
943computeFaceProperties(const Intersection& intersection,
944 const int,
945 const int,
946 const int,
947 const int,
948 DimVector& faceCenterInside,
949 DimVector& faceCenterOutside,
950 DimVector& faceAreaNormal,
951 /*isCpGrid=*/std::false_type) const
952{
953 // default implementation for DUNE grids
954 const auto& geometry = intersection.geometry();
955 faceCenterInside = geometry.center();
956 faceCenterOutside = faceCenterInside;
957
958 faceAreaNormal = intersection.centerUnitOuterNormal();
959 faceAreaNormal *= geometry.volume();
960}
961
962template<class Grid, class GridView, class ElementMapper, class CartesianIndexMapper, class Scalar>
963template<class Intersection>
965computeFaceProperties(const Intersection& 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 /*isCpGrid=*/std::true_type) const
974{
975 int faceIdx = intersection.id();
976
977 if(grid_.maxLevel() == 0) {
978 faceCenterInside = grid_.faceCenterEcl(insideElemIdx, insideFaceIdx, intersection);
979 faceCenterOutside = grid_.faceCenterEcl(outsideElemIdx, outsideFaceIdx, intersection);
980 faceAreaNormal = grid_.faceAreaNormalEcl(faceIdx);
981 }
982 else {
983 if ((intersection.inside().level() != intersection.outside().level())) {
984
985 // For CpGrid with LGRs, intersection laying on the boundary of an LGR, having two neighboring cells:
986 // one coarse neighboring cell and one refined neighboring cell, we get the corresponding parent
987 // intersection (from level 0), and use the center of the parent intersection for the coarse
988 // neighboring cell.
989
990 // Get parent intersection and its geometry
991 const auto& parentIntersection = grid_.getParentIntersectionFromLgrBoundaryFace(intersection);
992 const auto& parentIntersectionGeometry = parentIntersection.geometry();
993
994 // For the coarse neighboring cell, take the center of the parent intersection.
995 // For the refined neighboring cell, take the 'usual' center.
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);
1000
1001 // For some computations, it seems to be benefitial to replace the actual area of the refined face, by
1002 // the area of its parent face.
1003 // faceAreaNormal = parentIntersection.centerUnitOuterNormal();
1004 // faceAreaNormal *= parentIntersectionGeometry.volume();
1005
1007 faceAreaNormal = intersection.centerUnitOuterNormal();
1008 faceAreaNormal *= intersection.geometry().volume();
1009 }
1010 else {
1011 assert(intersection.inside().level() == intersection.outside().level());
1012
1013 faceCenterInside = grid_.faceCenterEcl(insideElemIdx, insideFaceIdx, intersection);
1014 faceCenterOutside = grid_.faceCenterEcl(outsideElemIdx, outsideFaceIdx, intersection);
1015
1016 // When the CpGrid has LGRs, we compute the face area normal differently.
1017 if (intersection.inside().level() > 0) { // remove intersection.inside().level() > 0
1018 faceAreaNormal = intersection.centerUnitOuterNormal();
1019 faceAreaNormal *= intersection.geometry().volume();
1020 }
1021 else {
1022 faceAreaNormal = grid_.faceAreaNormalEcl(faceIdx);
1023 }
1024 }
1025 }
1026}
1027template<class Grid, class GridView, class ElementMapper, class CartesianIndexMapper, class Scalar>
1028void
1030applyPinchNncToGridTrans_(const std::unordered_map<std::size_t,int>& cartesianToCompressed)
1031{
1032 // First scale NNCs with EDITNNC.
1033 const auto& nnc_input = eclState_.getPinchNNC();
1034
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;
1042
1043 if (low > high)
1044 std::swap(low, high);
1045
1046 if (low == -1 && high == -1)
1047 // Silently discard as it is not between active cells
1048 continue;
1049
1050 if (low == -1 || high == -1) {
1051 // We can end up here if one of the cells is overlap/ghost, because those
1052 // are lacking connections to other cells in the ghost/overlap.
1053 // Hence discard the NNC if it is between active cell and inactive cell
1054 continue;
1055 }
1056
1057 {
1058 auto candidate = trans_.find(details::isId(low, high));
1059 if (candidate != trans_.end()) {
1060 // the correctly calculated transmissibility is stored in
1061 // the NNC. Overwrite previous value with it.
1062 candidate->second = nncEntry.trans;
1063 }
1064 }
1065 }
1066}
1067
1068template<class Grid, class GridView, class ElementMapper, class CartesianIndexMapper, class Scalar>
1069void
1071applyNncToGridTrans_(const std::unordered_map<std::size_t,int>& cartesianToCompressed)
1072{
1073 // First scale NNCs with EDITNNC.
1074 const auto& nnc_input = eclState_.getInputNNC().input();
1075
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;
1083
1084 if (low > high)
1085 std::swap(low, high);
1086
1087 if (low == -1 && high == -1)
1088 // Silently discard as it is not between active cells
1089 continue;
1090
1091 if (low == -1 || high == -1) {
1092 // Discard the NNC if it is between active cell and inactive cell
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());
1097 continue;
1098 }
1099
1100 {
1101 auto candidate = trans_.find(details::isId(low, high));
1102 if (candidate != trans_.end()) {
1103 // NNC is represented by the grid and might be a neighboring connection
1104 // In this case the transmissibilty is added to the value already
1105 // set or computed.
1106 candidate->second += nncEntry.trans;
1107 }
1108 }
1109 // if (enableEnergy_) {
1110 // auto candidate = thermalHalfTrans_.find(details::directionalIsId(low, high));
1111 // if (candidate != trans_.end()) {
1112 // // NNC is represented by the grid and might be a neighboring connection
1113 // // In this case the transmissibilty is added to the value already
1114 // // set or computed.
1115 // candidate->second += nncEntry.transEnergy1;
1116 // }
1117 // auto candidate = thermalHalfTrans_.find(details::directionalIsId(high, low));
1118 // if (candidate != trans_.end()) {
1119 // // NNC is represented by the grid and might be a neighboring connection
1120 // // In this case the transmissibilty is added to the value already
1121 // // set or computed.
1122 // candidate->second += nncEntry.transEnergy2;
1123 // }
1124 // }
1125 // if (enableDiffusivity_) {
1126 // auto candidate = diffusivity_.find(details::isId(low, high));
1127 // if (candidate != trans_.end()) {
1128 // // NNC is represented by the grid and might be a neighboring connection
1129 // // In this case the transmissibilty is added to the value already
1130 // // set or computed.
1131 // candidate->second += nncEntry.transDiffusion;
1132 // }
1133 // }
1134 }
1135}
1136
1137template<class Grid, class GridView, class ElementMapper, class CartesianIndexMapper, class Scalar>
1139applyEditNncToGridTrans_(const std::unordered_map<std::size_t,int>& globalToLocal)
1140{
1141 const auto& input = eclState_.getInputNNC();
1142 applyEditNncToGridTransHelper_(globalToLocal, "EDITNNC",
1143 input.edit(),
1144 [&input](const NNCdata& nnc){
1145 return input.edit_location(nnc);},
1146 // Multiply transmissibility with EDITNNC value
1147 [](Scalar& trans, const Scalar& rhs){ trans *= rhs;});
1148}
1149
1150template<class Grid, class GridView, class ElementMapper, class CartesianIndexMapper, class Scalar>
1152applyEditNncrToGridTrans_(const std::unordered_map<std::size_t,int>& globalToLocal)
1153{
1154 const auto& input = eclState_.getInputNNC();
1155 applyEditNncToGridTransHelper_(globalToLocal, "EDITNNCR",
1156 input.editr(),
1157 [&input](const NNCdata& nnc){
1158 return input.editr_location(nnc);},
1159 // Replace Transmissibility with EDITNNCR value
1160 [](Scalar& trans, const Scalar& rhs){ trans = rhs;});
1161}
1162
1163template<class Grid, class GridView, class ElementMapper, class CartesianIndexMapper, class Scalar>
1165applyEditNncToGridTransHelper_(const std::unordered_map<std::size_t,int>& globalToLocal,
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)
1170{
1171 if (nncs.empty())
1172 return;
1173 const auto& cartDims = cartMapper_.cartesianDimensions();
1174
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];
1179
1180 return fmt::format("({},{},{})", i + 1,j + 1,k + 1);
1181 };
1182
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"
1186 "In {} line {} \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);
1190 };
1191
1192 // editNnc is supposed to only reference non-neighboring connections and not
1193 // neighboring connections. Use all entries for scaling if there is an NNC.
1194 // variable nnc incremented in loop body.
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);
1203
1204 if (lowIt == globalToLocal.end() || highIt == globalToLocal.end()) {
1205 // Prevent warnings for NNCs stored on other processes in parallel (both cells inactive)
1206 if ( lowIt != highIt && warnEditNNC_) {
1207 print_warning(*nnc);
1208 warning_count++;
1209 }
1210 ++nnc;
1211 continue;
1212 }
1213
1214 auto low = lowIt->second, high = highIt->second;
1215
1216 if (low > high)
1217 std::swap(low, high);
1218
1219 auto candidate = trans_.find(details::isId(low, high));
1220 if (candidate == trans_.end() && warnEditNNC_) {
1221 print_warning(*nnc);
1222 ++nnc;
1223 warning_count++;
1224 }
1225 else {
1226 // NNC exists
1227 while (nnc!= end && c1==nnc->cell1 && c2==nnc->cell2) {
1228 apply(candidate->second, nnc->trans);
1229 ++nnc;
1230 }
1231 }
1232 }
1233
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);
1238 }
1239}
1240
1241template<class Grid, class GridView, class ElementMapper, class CartesianIndexMapper, class Scalar>
1242void
1243Transmissibility<Grid,GridView,ElementMapper,CartesianIndexMapper,Scalar>::
1244applyNncMultreg_(const std::unordered_map<std::size_t,int>& cartesianToCompressed)
1245{
1246 const auto& inputNNC = this->eclState_.getInputNNC();
1247 const auto& transMult = this->eclState_.getTransMult();
1248
1249 auto compressedIdx = [&cartesianToCompressed](const std::size_t globIdx)
1250 {
1251 auto ixPos = cartesianToCompressed.find(globIdx);
1252 return (ixPos == cartesianToCompressed.end()) ? -1 : ixPos->second;
1253 };
1254
1255 // Apply region-based transmissibility multipliers (i.e., the MULTREGT
1256 // keyword) to those transmissibilities that are directly assigned from
1257 // the input.
1258 //
1259 // * NNC::input() covers the NNC keyword and any numerical aquifers
1260 // * NNC::editr() covers the EDITNNCR keyword
1261 //
1262 // Note: We do not apply MULTREGT to the entries in NNC::edit() since
1263 // those act as regular multipliers and have already been fully
1264 // accounted for in the multiplier part of the main loop of update() and
1265 // the applyEditNncToGridTrans_() member function.
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;
1270
1271 auto low = compressedIdx(c1);
1272 auto high = compressedIdx(c2);
1273
1274 if ((low == -1) || (high == -1)) {
1275 continue;
1276 }
1277
1278 if (low > high) {
1279 std::swap(low, high);
1280 }
1281
1282 auto candidate = this->trans_.find(details::isId(low, high));
1283 if (candidate != this->trans_.end()) {
1284 candidate->second *= transMult.getRegionMultiplierNNC(c1, c2);
1285 }
1286 }
1287 }
1288}
1289
1290template<class Grid, class GridView, class ElementMapper, class CartesianIndexMapper, class Scalar>
1291void Transmissibility<Grid,GridView,ElementMapper,CartesianIndexMapper,Scalar>::
1292computeHalfTrans_(Scalar& halfTrans,
1293 const DimVector& areaNormal,
1294 int faceIdx, // in the reference element that contains the intersection
1295 const DimVector& distance,
1296 const DimMatrix& perm) const
1297{
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();
1304}
1305
1306template<class Grid, class GridView, 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
1312{
1313 halfDiff = poro;
1314 halfDiff *= std::abs(Dune::dot(areaNormal, distance));
1315 halfDiff /= distance.two_norm2();
1316}
1317
1318template<class Grid, class GridView, 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
1323{
1324 const auto& cellCenter = centroids_(cellIdx);
1325 DimVector x = faceCenter;
1326 for (unsigned dimIdx = 0; dimIdx < dimWorld; ++dimIdx)
1327 x[dimIdx] -= cellCenter[dimIdx];
1328
1329 return x;
1330}
1331
1332template<class Grid, class GridView, class ElementMapper, class CartesianIndexMapper, class Scalar>
1333void Transmissibility<Grid,GridView,ElementMapper,CartesianIndexMapper,Scalar>::
1334applyMultipliers_(Scalar& trans,
1335 unsigned faceIdx,
1336 unsigned cartElemIdx,
1337 const TransMult& transMult) const
1338{
1339 // apply multiplyer for the transmissibility of the face. (the
1340 // face index is the index of the reference-element face which
1341 // contains the intersection of interest.)
1342 switch (faceIdx) {
1343 case 0: // left
1344 trans *= transMult.getMultiplier(cartElemIdx, FaceDir::XMinus);
1345 break;
1346 case 1: // right
1347 trans *= transMult.getMultiplier(cartElemIdx, FaceDir::XPlus);
1348 break;
1349
1350 case 2: // front
1351 trans *= transMult.getMultiplier(cartElemIdx, FaceDir::YMinus);
1352 break;
1353 case 3: // back
1354 trans *= transMult.getMultiplier(cartElemIdx, FaceDir::YPlus);
1355 break;
1356
1357 case 4: // bottom
1358 trans *= transMult.getMultiplier(cartElemIdx, FaceDir::ZMinus);
1359 break;
1360 case 5: // top
1361 trans *= transMult.getMultiplier(cartElemIdx, FaceDir::ZPlus);
1362 break;
1363 }
1364}
1365
1366template<class Grid, class GridView, class ElementMapper, class CartesianIndexMapper, class Scalar>
1367void Transmissibility<Grid,GridView,ElementMapper,CartesianIndexMapper,Scalar>::
1368applyNtg_(Scalar& trans,
1369 unsigned faceIdx,
1370 unsigned elemIdx,
1371 const std::vector<double>& ntg) const
1372{
1373 // apply multiplyer for the transmissibility of the face. (the
1374 // face index is the index of the reference-element face which
1375 // contains the intersection of interest.)
1376 switch (faceIdx) {
1377 case 0: // left
1378 trans *= ntg[elemIdx];
1379 break;
1380 case 1: // right
1381 trans *= ntg[elemIdx];
1382 break;
1383
1384 case 2: // front
1385 trans *= ntg[elemIdx];
1386 break;
1387 case 3: // back
1388 trans *= ntg[elemIdx];
1389 break;
1390
1391 // NTG does not apply to top and bottom faces
1392 }
1393}
1394
1395} // namespace Opm
1396
1397#endif // OPM_TRANSMISSIBILITY_IMPL_HPP
Definition Transmissibility.hpp:54
This file contains a set of helper functions used by VFPProd / VFPInj.
Definition blackoilboundaryratevector.hh:37
constexpr auto getPropValue()
get the value data member of a property
Definition propertysystem.hh:242