mirror of
https://github.com/OPM/opm-simulators.git
synced 2024-12-24 16:30:02 -06:00
1178 lines
49 KiB
C++
1178 lines
49 KiB
C++
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
|
|
// vi: set et ts=4 sw=4 sts=4:
|
|
/*
|
|
This file is part of the Open Porous Media project (OPM).
|
|
|
|
OPM is free software: you can redistribute it and/or modify
|
|
it under the terms of the GNU General Public License as published by
|
|
the Free Software Foundation, either version 2 of the License, or
|
|
(at your option) any later version.
|
|
|
|
OPM is distributed in the hope that it will be useful,
|
|
but WITHOUT ANY WARRANTY; without even the implied warranty of
|
|
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
|
GNU General Public License for more details.
|
|
|
|
You should have received a copy of the GNU General Public License
|
|
along with OPM. If not, see <http://www.gnu.org/licenses/>.
|
|
|
|
Consult the COPYING file in the top-level source directory of this
|
|
module for the precise wording of the license and the list of
|
|
copyright holders.
|
|
*/
|
|
|
|
#include <config.h>
|
|
#include <ebos/ecltransmissibility.hh>
|
|
|
|
#include <dune/common/version.hh>
|
|
#include <dune/grid/common/mcmgmapper.hh>
|
|
|
|
#include <opm/grid/CpGrid.hpp>
|
|
#include <opm/grid/polyhedralgrid.hh>
|
|
|
|
#if HAVE_DUNE_ALUGRID
|
|
#include <dune/alugrid/grid.hh>
|
|
#include <dune/alugrid/3d/gridview.hh>
|
|
#include "alucartesianindexmapper.hh"
|
|
#endif // HAVE_DUNE_ALUGRID
|
|
|
|
#include <opm/input/eclipse/EclipseState/EclipseState.hpp>
|
|
#include <opm/input/eclipse/EclipseState/Grid/FaceDir.hpp>
|
|
#include <opm/input/eclipse/EclipseState/Grid/FieldPropsManager.hpp>
|
|
#include <opm/input/eclipse/EclipseState/Grid/TransMult.hpp>
|
|
#include <opm/input/eclipse/Units/Units.hpp>
|
|
|
|
#if HAVE_DUNE_FEM
|
|
#include <dune/fem/gridpart/adaptiveleafgridpart.hh>
|
|
#include <dune/fem/gridpart/common/gridpart2gridview.hh>
|
|
#include <ebos/femcpgridcompat.hh>
|
|
#endif
|
|
|
|
#include <fmt/format.h>
|
|
|
|
#include <algorithm>
|
|
#include <cassert>
|
|
#include <cmath>
|
|
#include <sstream>
|
|
#include <stdexcept>
|
|
|
|
namespace {
|
|
|
|
constexpr unsigned elemIdxShift = 32; // bits
|
|
|
|
std::uint64_t isId(std::uint32_t elemIdx1, std::uint32_t elemIdx2)
|
|
{
|
|
std::uint32_t elemAIdx = std::min(elemIdx1, elemIdx2);
|
|
std::uint64_t elemBIdx = std::max(elemIdx1, elemIdx2);
|
|
|
|
return (elemBIdx<<elemIdxShift) + elemAIdx;
|
|
}
|
|
|
|
std::pair<std::uint32_t, std::uint32_t> isIdReverse(const std::uint64_t& id)
|
|
{
|
|
// Assigning an unsigned integer to a narrower type discards the most significant bits.
|
|
// See "The C programming language", section A.6.2.
|
|
// NOTE that the ordering of element A and B may have changed
|
|
std::uint32_t elemAIdx = id;
|
|
std::uint32_t elemBIdx = (id - elemAIdx) >> elemIdxShift;
|
|
|
|
return std::make_pair(elemAIdx, elemBIdx);
|
|
}
|
|
|
|
std::uint64_t directionalIsId(std::uint32_t elemIdx1, std::uint32_t elemIdx2)
|
|
{
|
|
return (std::uint64_t(elemIdx1)<<elemIdxShift) + elemIdx2;
|
|
}
|
|
|
|
}
|
|
|
|
namespace Opm {
|
|
|
|
template<class Grid, class GridView, class ElementMapper, class CartesianIndexMapper, class Scalar>
|
|
EclTransmissibility<Grid,GridView,ElementMapper,CartesianIndexMapper,Scalar>::
|
|
EclTransmissibility(const EclipseState& eclState,
|
|
const GridView& gridView,
|
|
const CartesianIndexMapper& cartMapper,
|
|
const Grid& grid,
|
|
std::function<std::array<double,dimWorld>(int)> centroids,
|
|
bool enableEnergy,
|
|
bool enableDiffusivity)
|
|
: eclState_(eclState)
|
|
, gridView_(gridView)
|
|
, cartMapper_(cartMapper)
|
|
, grid_(grid)
|
|
, centroids_(centroids)
|
|
, enableEnergy_(enableEnergy)
|
|
, enableDiffusivity_(enableDiffusivity)
|
|
{
|
|
const UnitSystem& unitSystem = eclState_.getDeckUnitSystem();
|
|
transmissibilityThreshold_ = unitSystem.parse("Transmissibility").getSIScaling() * 1e-6;
|
|
}
|
|
|
|
template<class Grid, class GridView, class ElementMapper, class CartesianIndexMapper, class Scalar>
|
|
Scalar EclTransmissibility<Grid,GridView,ElementMapper,CartesianIndexMapper,Scalar>::
|
|
transmissibility(unsigned elemIdx1, unsigned elemIdx2) const
|
|
{
|
|
return trans_.at(isId(elemIdx1, elemIdx2));
|
|
}
|
|
|
|
template<class Grid, class GridView, class ElementMapper, class CartesianIndexMapper, class Scalar>
|
|
Scalar EclTransmissibility<Grid,GridView,ElementMapper,CartesianIndexMapper,Scalar>::
|
|
transmissibilityBoundary(unsigned elemIdx, unsigned boundaryFaceIdx) const
|
|
{
|
|
return transBoundary_.at(std::make_pair(elemIdx, boundaryFaceIdx));
|
|
}
|
|
|
|
template<class Grid, class GridView, class ElementMapper, class CartesianIndexMapper, class Scalar>
|
|
Scalar EclTransmissibility<Grid,GridView,ElementMapper,CartesianIndexMapper,Scalar>::
|
|
thermalHalfTrans(unsigned insideElemIdx, unsigned outsideElemIdx) const
|
|
{
|
|
return thermalHalfTrans_.at(directionalIsId(insideElemIdx, outsideElemIdx));
|
|
}
|
|
|
|
template<class Grid, class GridView, class ElementMapper, class CartesianIndexMapper, class Scalar>
|
|
Scalar EclTransmissibility<Grid,GridView,ElementMapper,CartesianIndexMapper,Scalar>::
|
|
thermalHalfTransBoundary(unsigned insideElemIdx, unsigned boundaryFaceIdx) const
|
|
{
|
|
return thermalHalfTransBoundary_.at(std::make_pair(insideElemIdx, boundaryFaceIdx));
|
|
}
|
|
|
|
template<class Grid, class GridView, class ElementMapper, class CartesianIndexMapper, class Scalar>
|
|
Scalar EclTransmissibility<Grid,GridView,ElementMapper,CartesianIndexMapper,Scalar>::
|
|
diffusivity(unsigned elemIdx1, unsigned elemIdx2) const
|
|
{
|
|
if (diffusivity_.empty())
|
|
return 0.0;
|
|
|
|
return diffusivity_.at(isId(elemIdx1, elemIdx2));
|
|
|
|
}
|
|
|
|
template<class Grid, class GridView, class ElementMapper, class CartesianIndexMapper, class Scalar>
|
|
void EclTransmissibility<Grid,GridView,ElementMapper,CartesianIndexMapper,Scalar>::
|
|
update(bool global, const std::function<unsigned int(unsigned int)>& map)
|
|
{
|
|
const auto& cartDims = cartMapper_.cartesianDimensions();
|
|
auto& transMult = eclState_.getTransMult();
|
|
const auto& comm = gridView_.comm();
|
|
ElementMapper elemMapper(gridView_, Dune::mcmgElementLayout());
|
|
|
|
// get the ntg values, the ntg values are modified for the cells merged with minpv
|
|
const std::vector<double>& ntg = eclState_.fieldProps().get_double("NTG");
|
|
const bool updateDiffusivity = eclState_.getSimulationConfig().isDiffusive() && enableDiffusivity_;
|
|
unsigned numElements = elemMapper.size();
|
|
|
|
if (map)
|
|
extractPermeability_(map);
|
|
else
|
|
extractPermeability_();
|
|
|
|
// calculate the axis specific centroids of all elements
|
|
std::array<std::vector<DimVector>, dimWorld> axisCentroids;
|
|
|
|
for (unsigned dimIdx = 0; dimIdx < dimWorld; ++dimIdx)
|
|
axisCentroids[dimIdx].resize(numElements);
|
|
|
|
auto elemIt = gridView_.template begin</*codim=*/ 0>();
|
|
const auto& elemEndIt = gridView_.template end</*codim=*/ 0>();
|
|
size_t centroidIdx = 0;
|
|
for (; elemIt != elemEndIt; ++elemIt, ++centroidIdx) {
|
|
const auto& elem = *elemIt;
|
|
unsigned elemIdx = elemMapper.index(elem);
|
|
|
|
// compute the axis specific "centroids" used for the transmissibilities. for
|
|
// consistency with the flow simulator, we use the element centers as
|
|
// computed by opm-parser's Opm::EclipseGrid class for all axes.
|
|
std::array<double, dimWorld> centroid = centroids_(elemIdx);
|
|
|
|
for (unsigned axisIdx = 0; axisIdx < dimWorld; ++axisIdx)
|
|
for (unsigned dimIdx = 0; dimIdx < dimWorld; ++dimIdx)
|
|
axisCentroids[axisIdx][elemIdx][dimIdx] = centroid[dimIdx];
|
|
}
|
|
|
|
// reserving some space in the hashmap upfront saves quite a bit of time because
|
|
// resizes are costly for hashmaps and there would be quite a few of them if we
|
|
// would not have a rough idea of how large the final map will be (the rough idea
|
|
// is a conforming Cartesian grid).
|
|
trans_.clear();
|
|
trans_.reserve(numElements*3*1.05);
|
|
|
|
transBoundary_.clear();
|
|
|
|
// if energy is enabled, let's do the same for the "thermal half transmissibilities"
|
|
if (enableEnergy_) {
|
|
thermalHalfTrans_.clear();
|
|
thermalHalfTrans_.reserve(numElements*6*1.05);
|
|
|
|
thermalHalfTransBoundary_.clear();
|
|
}
|
|
|
|
// if diffusion is enabled, let's do the same for the "diffusivity"
|
|
if (updateDiffusivity) {
|
|
diffusivity_.clear();
|
|
diffusivity_.reserve(numElements*3*1.05);
|
|
extractPorosity_();
|
|
}
|
|
|
|
// The MULTZ needs special case if the option is ALL
|
|
// Then the smallest multiplier is applied.
|
|
// Default is to apply the top and bottom multiplier
|
|
bool useSmallestMultiplier;
|
|
if (comm.rank() == 0) {
|
|
const auto& eclGrid = eclState_.getInputGrid();
|
|
useSmallestMultiplier = eclGrid.getMultzOption() == PinchMode::ModeEnum::ALL;
|
|
}
|
|
if (global && comm.size() > 1) {
|
|
comm.broadcast(&useSmallestMultiplier, 1, 0);
|
|
}
|
|
|
|
// compute the transmissibilities for all intersections
|
|
elemIt = gridView_.template begin</*codim=*/ 0>();
|
|
for (; elemIt != elemEndIt; ++elemIt) {
|
|
const auto& elem = *elemIt;
|
|
unsigned elemIdx = elemMapper.index(elem);
|
|
|
|
auto isIt = gridView_.ibegin(elem);
|
|
const auto& isEndIt = gridView_.iend(elem);
|
|
unsigned boundaryIsIdx = 0;
|
|
for (; isIt != isEndIt; ++ isIt) {
|
|
// store intersection, this might be costly
|
|
const auto& intersection = *isIt;
|
|
|
|
// deal with grid boundaries
|
|
if (intersection.boundary()) {
|
|
// compute the transmissibilty for the boundary intersection
|
|
const auto& geometry = intersection.geometry();
|
|
const auto& faceCenterInside = geometry.center();
|
|
|
|
auto faceAreaNormal = intersection.centerUnitOuterNormal();
|
|
faceAreaNormal *= geometry.volume();
|
|
|
|
Scalar transBoundaryIs;
|
|
computeHalfTrans_(transBoundaryIs,
|
|
faceAreaNormal,
|
|
intersection.indexInInside(),
|
|
distanceVector_(faceCenterInside,
|
|
intersection.indexInInside(),
|
|
elemIdx,
|
|
axisCentroids),
|
|
permeability_[elemIdx]);
|
|
|
|
// normally there would be two half-transmissibilities that would be
|
|
// averaged. on the grid boundary there only is the half
|
|
// transmissibility of the interior element.
|
|
transBoundary_[std::make_pair(elemIdx, boundaryIsIdx)] = transBoundaryIs;
|
|
|
|
// for boundary intersections we also need to compute the thermal
|
|
// half transmissibilities
|
|
if (enableEnergy_) {
|
|
Scalar transBoundaryEnergyIs;
|
|
computeHalfDiffusivity_(transBoundaryEnergyIs,
|
|
faceAreaNormal,
|
|
distanceVector_(faceCenterInside,
|
|
intersection.indexInInside(),
|
|
elemIdx,
|
|
axisCentroids),
|
|
1.0);
|
|
thermalHalfTransBoundary_[std::make_pair(elemIdx, boundaryIsIdx)] =
|
|
transBoundaryEnergyIs;
|
|
}
|
|
|
|
++ boundaryIsIdx;
|
|
continue;
|
|
}
|
|
|
|
if (!intersection.neighbor()) {
|
|
// elements can be on process boundaries, i.e. they are not on the
|
|
// domain boundary yet they don't have neighbors.
|
|
++ boundaryIsIdx;
|
|
continue;
|
|
}
|
|
|
|
const auto& outsideElem = intersection.outside();
|
|
unsigned outsideElemIdx = elemMapper.index(outsideElem);
|
|
|
|
unsigned insideCartElemIdx = cartMapper_.cartesianIndex(elemIdx);
|
|
unsigned outsideCartElemIdx = cartMapper_.cartesianIndex(outsideElemIdx);
|
|
|
|
// we only need to calculate a face's transmissibility
|
|
// once...
|
|
if (insideCartElemIdx > outsideCartElemIdx)
|
|
continue;
|
|
|
|
// local indices of the faces of the inside and
|
|
// outside elements which contain the intersection
|
|
int insideFaceIdx = intersection.indexInInside();
|
|
int outsideFaceIdx = intersection.indexInOutside();
|
|
|
|
if (insideFaceIdx == -1) {
|
|
// NNC. Set zero transmissibility, as it will be
|
|
// *added to* by applyNncToGridTrans_() later.
|
|
assert(outsideFaceIdx == -1);
|
|
trans_[isId(elemIdx, outsideElemIdx)] = 0.0;
|
|
continue;
|
|
}
|
|
|
|
DimVector faceCenterInside;
|
|
DimVector faceCenterOutside;
|
|
DimVector faceAreaNormal;
|
|
|
|
typename std::is_same<Grid, Dune::CpGrid>::type isCpGrid;
|
|
computeFaceProperties(intersection,
|
|
elemIdx,
|
|
insideFaceIdx,
|
|
outsideElemIdx,
|
|
outsideFaceIdx,
|
|
faceCenterInside,
|
|
faceCenterOutside,
|
|
faceAreaNormal,
|
|
isCpGrid);
|
|
|
|
Scalar halfTrans1;
|
|
Scalar halfTrans2;
|
|
|
|
computeHalfTrans_(halfTrans1,
|
|
faceAreaNormal,
|
|
insideFaceIdx,
|
|
distanceVector_(faceCenterInside,
|
|
intersection.indexInInside(),
|
|
elemIdx,
|
|
axisCentroids),
|
|
permeability_[elemIdx]);
|
|
computeHalfTrans_(halfTrans2,
|
|
faceAreaNormal,
|
|
outsideFaceIdx,
|
|
distanceVector_(faceCenterOutside,
|
|
intersection.indexInOutside(),
|
|
outsideElemIdx,
|
|
axisCentroids),
|
|
permeability_[outsideElemIdx]);
|
|
|
|
applyNtg_(halfTrans1, insideFaceIdx, elemIdx, ntg);
|
|
applyNtg_(halfTrans2, outsideFaceIdx, outsideElemIdx, ntg);
|
|
|
|
// convert half transmissibilities to full face
|
|
// transmissibilities using the harmonic mean
|
|
Scalar trans;
|
|
if (std::abs(halfTrans1) < 1e-30 || std::abs(halfTrans2) < 1e-30)
|
|
// avoid division by zero
|
|
trans = 0.0;
|
|
else
|
|
trans = 1.0 / (1.0/halfTrans1 + 1.0/halfTrans2);
|
|
|
|
// apply the full face transmissibility multipliers
|
|
// for the inside ...
|
|
|
|
if (useSmallestMultiplier)
|
|
{
|
|
// Currently PINCH(4) is never queries and hence PINCH(4) == TOPBOT is assumed
|
|
// and in this branch PINCH(5) == ALL holds
|
|
applyAllZMultipliers_(trans, insideFaceIdx, outsideFaceIdx, insideCartElemIdx,
|
|
outsideCartElemIdx, transMult, cartDims,
|
|
/* pinchTop= */ false);
|
|
}
|
|
else
|
|
{
|
|
applyMultipliers_(trans, insideFaceIdx, insideCartElemIdx, transMult);
|
|
// ... and outside elements
|
|
applyMultipliers_(trans, outsideFaceIdx, outsideCartElemIdx, transMult);
|
|
}
|
|
|
|
// apply the region multipliers (cf. the MULTREGT keyword)
|
|
FaceDir::DirEnum faceDir;
|
|
switch (insideFaceIdx) {
|
|
case 0:
|
|
case 1:
|
|
faceDir = FaceDir::XPlus;
|
|
break;
|
|
|
|
case 2:
|
|
case 3:
|
|
faceDir = FaceDir::YPlus;
|
|
break;
|
|
|
|
case 4:
|
|
case 5:
|
|
faceDir = FaceDir::ZPlus;
|
|
break;
|
|
|
|
default:
|
|
throw std::logic_error("Could not determine a face direction");
|
|
}
|
|
|
|
trans *= transMult.getRegionMultiplier(insideCartElemIdx,
|
|
outsideCartElemIdx,
|
|
faceDir);
|
|
|
|
trans_[isId(elemIdx, outsideElemIdx)] = trans;
|
|
|
|
// update the "thermal half transmissibility" for the intersection
|
|
if (enableEnergy_) {
|
|
|
|
Scalar halfDiffusivity1;
|
|
Scalar halfDiffusivity2;
|
|
|
|
computeHalfDiffusivity_(halfDiffusivity1,
|
|
faceAreaNormal,
|
|
distanceVector_(faceCenterInside,
|
|
intersection.indexInInside(),
|
|
elemIdx,
|
|
axisCentroids),
|
|
1.0);
|
|
computeHalfDiffusivity_(halfDiffusivity2,
|
|
faceAreaNormal,
|
|
distanceVector_(faceCenterOutside,
|
|
intersection.indexInOutside(),
|
|
outsideElemIdx,
|
|
axisCentroids),
|
|
1.0);
|
|
//TODO Add support for multipliers
|
|
thermalHalfTrans_[directionalIsId(elemIdx, outsideElemIdx)] = halfDiffusivity1;
|
|
thermalHalfTrans_[directionalIsId(outsideElemIdx, elemIdx)] = halfDiffusivity2;
|
|
}
|
|
|
|
// update the "diffusive half transmissibility" for the intersection
|
|
if (updateDiffusivity) {
|
|
|
|
Scalar halfDiffusivity1;
|
|
Scalar halfDiffusivity2;
|
|
|
|
computeHalfDiffusivity_(halfDiffusivity1,
|
|
faceAreaNormal,
|
|
distanceVector_(faceCenterInside,
|
|
intersection.indexInInside(),
|
|
elemIdx,
|
|
axisCentroids),
|
|
porosity_[elemIdx]);
|
|
computeHalfDiffusivity_(halfDiffusivity2,
|
|
faceAreaNormal,
|
|
distanceVector_(faceCenterOutside,
|
|
intersection.indexInOutside(),
|
|
outsideElemIdx,
|
|
axisCentroids),
|
|
porosity_[outsideElemIdx]);
|
|
|
|
applyNtg_(halfDiffusivity1, insideFaceIdx, elemIdx, ntg);
|
|
applyNtg_(halfDiffusivity2, outsideFaceIdx, outsideElemIdx, ntg);
|
|
|
|
//TODO Add support for multipliers
|
|
Scalar diffusivity;
|
|
if (std::abs(halfDiffusivity1) < 1e-30 || std::abs(halfDiffusivity2) < 1e-30)
|
|
// avoid division by zero
|
|
diffusivity = 0.0;
|
|
else
|
|
diffusivity = 1.0 / (1.0/halfDiffusivity1 + 1.0/halfDiffusivity2);
|
|
|
|
|
|
diffusivity_[isId(elemIdx, outsideElemIdx)] = diffusivity;
|
|
}
|
|
}
|
|
}
|
|
|
|
// potentially overwrite and/or modify transmissibilities based on input from deck
|
|
updateFromEclState_(global);
|
|
|
|
// Create mapping from global to local index
|
|
std::unordered_map<std::size_t,int> globalToLocal;
|
|
|
|
// loop over all elements (global grid) and store Cartesian index
|
|
elemIt = grid_.leafGridView().template begin<0>();
|
|
|
|
for (; elemIt != elemEndIt; ++elemIt) {
|
|
int elemIdx = elemMapper.index(*elemIt);
|
|
int cartElemIdx = cartMapper_.cartesianIndex(elemIdx);
|
|
globalToLocal[cartElemIdx] = elemIdx;
|
|
}
|
|
applyEditNncToGridTrans_(globalToLocal);
|
|
applyNncToGridTrans_(globalToLocal);
|
|
|
|
//remove very small non-neighbouring transmissibilities
|
|
removeSmallNonCartesianTransmissibilities_();
|
|
}
|
|
|
|
template<class Grid, class GridView, class ElementMapper, class CartesianIndexMapper, class Scalar>
|
|
void EclTransmissibility<Grid,GridView,ElementMapper,CartesianIndexMapper,Scalar>::
|
|
extractPermeability_()
|
|
{
|
|
unsigned numElem = gridView_.size(/*codim=*/0);
|
|
permeability_.resize(numElem);
|
|
|
|
// read the intrinsic permeabilities from the eclState. Note that all arrays
|
|
// provided by eclState are one-per-cell of "uncompressed" grid, whereas the
|
|
// simulation grid might remove a few elements. (e.g. because it is distributed
|
|
// over several processes.)
|
|
const auto& fp = eclState_.fieldProps();
|
|
if (fp.has_double("PERMX")) {
|
|
const std::vector<double>& permxData = fp.get_double("PERMX");
|
|
|
|
std::vector<double> permyData;
|
|
if (fp.has_double("PERMY"))
|
|
permyData = fp.get_double("PERMY");
|
|
else
|
|
permyData = permxData;
|
|
|
|
std::vector<double> permzData;
|
|
if (fp.has_double("PERMZ"))
|
|
permzData = fp.get_double("PERMZ");
|
|
else
|
|
permzData = permxData;
|
|
|
|
for (size_t dofIdx = 0; dofIdx < numElem; ++ dofIdx) {
|
|
permeability_[dofIdx] = 0.0;
|
|
permeability_[dofIdx][0][0] = permxData[dofIdx];
|
|
permeability_[dofIdx][1][1] = permyData[dofIdx];
|
|
permeability_[dofIdx][2][2] = permzData[dofIdx];
|
|
}
|
|
|
|
// for now we don't care about non-diagonal entries
|
|
|
|
}
|
|
else
|
|
throw std::logic_error("Can't read the intrinsic permeability from the ecl state. "
|
|
"(The PERM{X,Y,Z} keywords are missing)");
|
|
}
|
|
|
|
template<class Grid, class GridView, class ElementMapper, class CartesianIndexMapper, class Scalar>
|
|
void EclTransmissibility<Grid,GridView,ElementMapper,CartesianIndexMapper,Scalar>::
|
|
extractPermeability_(const std::function<unsigned int(unsigned int)>& map)
|
|
{
|
|
unsigned numElem = gridView_.size(/*codim=*/0);
|
|
permeability_.resize(numElem);
|
|
|
|
// read the intrinsic permeabilities from the eclState. Note that all arrays
|
|
// provided by eclState are one-per-cell of "uncompressed" grid, whereas the
|
|
// simulation grid might remove a few elements. (e.g. because it is distributed
|
|
// over several processes.)
|
|
const auto& fp = eclState_.fieldProps();
|
|
if (fp.has_double("PERMX")) {
|
|
const std::vector<double>& permxData = fp.get_double("PERMX");
|
|
|
|
std::vector<double> permyData;
|
|
if (fp.has_double("PERMY"))
|
|
permyData = fp.get_double("PERMY");
|
|
else
|
|
permyData = permxData;
|
|
|
|
std::vector<double> permzData;
|
|
if (fp.has_double("PERMZ"))
|
|
permzData = fp.get_double("PERMZ");
|
|
else
|
|
permzData = permxData;
|
|
|
|
for (size_t dofIdx = 0; dofIdx < numElem; ++ dofIdx) {
|
|
permeability_[dofIdx] = 0.0;
|
|
size_t inputDofIdx = map(dofIdx);
|
|
permeability_[dofIdx][0][0] = permxData[inputDofIdx];
|
|
permeability_[dofIdx][1][1] = permyData[inputDofIdx];
|
|
permeability_[dofIdx][2][2] = permzData[inputDofIdx];
|
|
}
|
|
|
|
// for now we don't care about non-diagonal entries
|
|
|
|
}
|
|
else
|
|
throw std::logic_error("Can't read the intrinsic permeability from the ecl state. "
|
|
"(The PERM{X,Y,Z} keywords are missing)");
|
|
}
|
|
|
|
template<class Grid, class GridView, class ElementMapper, class CartesianIndexMapper, class Scalar>
|
|
void EclTransmissibility<Grid,GridView,ElementMapper,CartesianIndexMapper,Scalar>::
|
|
extractPorosity_()
|
|
{
|
|
// read the intrinsic porosity from the eclState. Note that all arrays
|
|
// provided by eclState are one-per-cell of "uncompressed" grid, whereas the
|
|
// simulation grid might remove a few elements. (e.g. because it is distributed
|
|
// over several processes.)
|
|
const auto& fp = eclState_.fieldProps();
|
|
if (fp.has_double("PORO")) {
|
|
porosity_ = fp.get_double("PORO");
|
|
}
|
|
else
|
|
throw std::logic_error("Can't read the porosityfrom the ecl state. "
|
|
"(The PORO keywords are missing)");
|
|
}
|
|
|
|
template<class Grid, class GridView, class ElementMapper, class CartesianIndexMapper, class Scalar>
|
|
void EclTransmissibility<Grid,GridView,ElementMapper,CartesianIndexMapper,Scalar>::
|
|
removeSmallNonCartesianTransmissibilities_()
|
|
{
|
|
const auto& cartDims = cartMapper_.cartesianDimensions();
|
|
for (auto&& trans: trans_) {
|
|
if (trans.second < transmissibilityThreshold_) {
|
|
const auto& id = trans.first;
|
|
const auto& elements = isIdReverse(id);
|
|
int gc1 = std::min(cartMapper_.cartesianIndex(elements.first), cartMapper_.cartesianIndex(elements.second));
|
|
int gc2 = std::max(cartMapper_.cartesianIndex(elements.first), cartMapper_.cartesianIndex(elements.second));
|
|
|
|
// only adjust the NNCs
|
|
if (gc2 - gc1 == 1 || gc2 - gc1 == cartDims[0] || gc2 - gc1 == cartDims[0]*cartDims[1])
|
|
continue;
|
|
|
|
//remove transmissibilities less than the threshold (by default 1e-6 in the deck's unit system)
|
|
trans.second = 0.0;
|
|
}
|
|
}
|
|
}
|
|
|
|
template<class Grid, class GridView, class ElementMapper, class CartesianIndexMapper, class Scalar>
|
|
void EclTransmissibility<Grid,GridView,ElementMapper,CartesianIndexMapper, Scalar>::
|
|
applyAllZMultipliers_(Scalar& trans,
|
|
unsigned insideFaceIdx,
|
|
unsigned outsideFaceIdx,
|
|
unsigned insideCartElemIdx,
|
|
unsigned outsideCartElemIdx,
|
|
const TransMult& transMult,
|
|
const std::array<int, dimWorld>& cartDims,
|
|
bool pinchTop)
|
|
{
|
|
if (insideFaceIdx > 3) { // top or or bottom
|
|
assert(insideFaceIdx==5); // as insideCartElemIdx < outsideCartElemIdx holds for the Z column
|
|
assert(outsideCartElemIdx > insideCartElemIdx);
|
|
auto lastCartElemIdx = outsideCartElemIdx - cartDims[0]*cartDims[1];
|
|
// Last multiplier
|
|
Scalar mult = transMult.getMultiplier(lastCartElemIdx , FaceDir::ZPlus);
|
|
|
|
if ( !pinchTop )
|
|
{
|
|
// pick the smallest multiplier for Z+ while looking down the pillar until reaching the other end of the connection
|
|
// While Z- is not all used here.
|
|
for(auto cartElemIdx = insideCartElemIdx; cartElemIdx < lastCartElemIdx;
|
|
cartElemIdx += cartDims[0]*cartDims[1])
|
|
{
|
|
mult = std::min(mult, transMult.getMultiplier(cartElemIdx, FaceDir::ZPlus));
|
|
}
|
|
}
|
|
|
|
trans *= mult;
|
|
applyMultipliers_(trans, outsideFaceIdx, outsideCartElemIdx, transMult);
|
|
}
|
|
else
|
|
{
|
|
applyMultipliers_(trans, insideFaceIdx, insideCartElemIdx, transMult);
|
|
applyMultipliers_(trans, outsideFaceIdx, outsideCartElemIdx, transMult);
|
|
}
|
|
}
|
|
|
|
template<class Grid, class GridView, class ElementMapper, class CartesianIndexMapper, class Scalar>
|
|
void EclTransmissibility<Grid,GridView,ElementMapper,CartesianIndexMapper,Scalar>::
|
|
updateFromEclState_(bool global)
|
|
{
|
|
const FieldPropsManager* fp =
|
|
(global) ? &(eclState_.fieldProps()) :
|
|
&(eclState_.globalFieldProps());
|
|
|
|
std::array<bool,3> is_tran {fp->tran_active("TRANX"),
|
|
fp->tran_active("TRANY"),
|
|
fp->tran_active("TRANZ")};
|
|
|
|
if( !(is_tran[0] ||is_tran[1] || is_tran[2]) )
|
|
{
|
|
// Skip unneeded expensive traversals
|
|
return;
|
|
}
|
|
|
|
std::array<std::string, 3> keywords {"TRANX", "TRANY", "TRANZ"};
|
|
std::array<std::vector<double>,3> trans = createTransmissibilityArrays_(is_tran);
|
|
auto key = keywords.begin();
|
|
auto perform = is_tran.begin();
|
|
|
|
for (auto it = trans.begin(); it != trans.end(); ++it, ++key, ++perform)
|
|
{
|
|
if(perform)
|
|
fp->apply_tran(*key, *it);
|
|
}
|
|
|
|
resetTransmissibilityFromArrays_(is_tran, trans);
|
|
}
|
|
|
|
template<class Grid, class GridView, class ElementMapper, class CartesianIndexMapper, class Scalar>
|
|
std::array<std::vector<double>,3>
|
|
EclTransmissibility<Grid,GridView,ElementMapper,CartesianIndexMapper,Scalar>::
|
|
createTransmissibilityArrays_(const std::array<bool,3>& is_tran)
|
|
{
|
|
const auto& cartDims = cartMapper_.cartesianDimensions();
|
|
ElementMapper elemMapper(gridView_, Dune::mcmgElementLayout());
|
|
|
|
auto numElem = gridView_.size(/*codim=*/0);
|
|
std::array<std::vector<double>,3> trans =
|
|
{ std::vector<double>(is_tran[0] ? numElem : 0, 0),
|
|
std::vector<double>(is_tran[1] ? numElem : 0, 0),
|
|
std::vector<double>(is_tran[2] ? numElem : 0, 0)};
|
|
|
|
// compute the transmissibilities for all intersections
|
|
auto elemIt = gridView_.template begin</*codim=*/ 0>();
|
|
const auto& elemEndIt = gridView_.template end</*codim=*/ 0>();
|
|
|
|
for (; elemIt != elemEndIt; ++elemIt) {
|
|
const auto& elem = *elemIt;
|
|
auto isIt = gridView_.ibegin(elem);
|
|
const auto& isEndIt = gridView_.iend(elem);
|
|
for (; isIt != isEndIt; ++ isIt) {
|
|
// store intersection, this might be costly
|
|
const auto& intersection = *isIt;
|
|
if (!intersection.neighbor())
|
|
continue; // intersection is on the domain boundary
|
|
|
|
// In the EclState TRANX[c1] is transmissibility in X+
|
|
// direction. Ordering of compressed (c1,c2) and cartesian index
|
|
// (gc1, gc2) is coherent (c1 < c2 <=> gc1 < gc2). This also
|
|
// holds for the global grid. While distributing changes the
|
|
// order of the local indices, the transmissibilities are still
|
|
// stored at the cell with the lower global cartesian index as
|
|
// the fieldprops are communicated by the grid.
|
|
unsigned c1 = elemMapper.index(intersection.inside());
|
|
unsigned c2 = elemMapper.index(intersection.outside());
|
|
int gc1 = cartMapper_.cartesianIndex(c1);
|
|
int gc2 = cartMapper_.cartesianIndex(c2);
|
|
|
|
if (gc1 > gc2)
|
|
continue; // we only need to handle each connection once, thank you.
|
|
|
|
auto isID = isId(c1, c2);
|
|
|
|
if (gc2 - gc1 == 1 && cartDims[0] > 1) {
|
|
if (is_tran[0])
|
|
// set simulator internal transmissibilities to values from inputTranx
|
|
trans[0][c1] = trans_[isID];
|
|
}
|
|
else if (gc2 - gc1 == cartDims[0] && cartDims[1] > 1) {
|
|
if (is_tran[1])
|
|
// set simulator internal transmissibilities to values from inputTrany
|
|
trans[1][c1] = trans_[isID];
|
|
}
|
|
else if (gc2 - gc1 == cartDims[0]*cartDims[1]) {
|
|
if (is_tran[2])
|
|
// set simulator internal transmissibilities to values from inputTranz
|
|
trans[2][c1] = trans_[isID];
|
|
}
|
|
//else.. We don't support modification of NNC at the moment.
|
|
}
|
|
}
|
|
|
|
return trans;
|
|
}
|
|
|
|
template<class Grid, class GridView, class ElementMapper, class CartesianIndexMapper, class Scalar>
|
|
void EclTransmissibility<Grid,GridView,ElementMapper,CartesianIndexMapper,Scalar>::
|
|
resetTransmissibilityFromArrays_(const std::array<bool,3>& is_tran,
|
|
const std::array<std::vector<double>,3>& trans)
|
|
{
|
|
const auto& cartDims = cartMapper_.cartesianDimensions();
|
|
ElementMapper elemMapper(gridView_, Dune::mcmgElementLayout());
|
|
|
|
// compute the transmissibilities for all intersections
|
|
auto elemIt = gridView_.template begin</*codim=*/ 0>();
|
|
const auto& elemEndIt = gridView_.template end</*codim=*/ 0>();
|
|
|
|
for (; elemIt != elemEndIt; ++elemIt) {
|
|
const auto& elem = *elemIt;
|
|
auto isIt = gridView_.ibegin(elem);
|
|
const auto& isEndIt = gridView_.iend(elem);
|
|
for (; isIt != isEndIt; ++ isIt) {
|
|
// store intersection, this might be costly
|
|
const auto& intersection = *isIt;
|
|
if (!intersection.neighbor())
|
|
continue; // intersection is on the domain boundary
|
|
|
|
// In the EclState TRANX[c1] is transmissibility in X+
|
|
// direction. Ordering of compressed (c1,c2) and cartesian index
|
|
// (gc1, gc2) is coherent (c1 < c2 <=> gc1 < gc2). This also
|
|
// holds for the global grid. While distributing changes the
|
|
// order of the local indices, the transmissibilities are still
|
|
// stored at the cell with the lower global cartesian index as
|
|
// the fieldprops are communicated by the grid.
|
|
unsigned c1 = elemMapper.index(intersection.inside());
|
|
unsigned c2 = elemMapper.index(intersection.outside());
|
|
int gc1 = cartMapper_.cartesianIndex(c1);
|
|
int gc2 = cartMapper_.cartesianIndex(c2);
|
|
if (gc1 > gc2)
|
|
continue; // we only need to handle each connection once, thank you.
|
|
|
|
auto isID = isId(c1, c2);
|
|
|
|
if (gc2 - gc1 == 1 && cartDims[0] > 1) {
|
|
if (is_tran[0])
|
|
// set simulator internal transmissibilities to values from inputTranx
|
|
trans_[isID] = trans[0][c1];
|
|
}
|
|
else if (gc2 - gc1 == cartDims[0] && cartDims[1] > 1) {
|
|
if (is_tran[1])
|
|
// set simulator internal transmissibilities to values from inputTrany
|
|
trans_[isID] = trans[1][c1];
|
|
}
|
|
else if (gc2 - gc1 == cartDims[0]*cartDims[1]) {
|
|
if (is_tran[2])
|
|
// set simulator internal transmissibilities to values from inputTranz
|
|
trans_[isID] = trans[2][c1];
|
|
}
|
|
//else.. We don't support modification of NNC at the moment.
|
|
}
|
|
}
|
|
}
|
|
|
|
template<class Grid, class GridView, class ElementMapper, class CartesianIndexMapper, class Scalar>
|
|
template<class Intersection>
|
|
void EclTransmissibility<Grid,GridView,ElementMapper,CartesianIndexMapper,Scalar>::
|
|
computeFaceProperties(const Intersection& intersection,
|
|
const int,
|
|
const int,
|
|
const int,
|
|
const int,
|
|
DimVector& faceCenterInside,
|
|
DimVector& faceCenterOutside,
|
|
DimVector& faceAreaNormal,
|
|
/*isCpGrid=*/std::false_type) const
|
|
{
|
|
// default implementation for DUNE grids
|
|
const auto& geometry = intersection.geometry();
|
|
faceCenterInside = geometry.center();
|
|
faceCenterOutside = faceCenterInside;
|
|
|
|
faceAreaNormal = intersection.centerUnitOuterNormal();
|
|
faceAreaNormal *= geometry.volume();
|
|
}
|
|
|
|
|
|
template<class Grid, class GridView, class ElementMapper, class CartesianIndexMapper, class Scalar>
|
|
template<class Intersection>
|
|
void EclTransmissibility<Grid,GridView,ElementMapper,CartesianIndexMapper,Scalar>::
|
|
computeFaceProperties(const Intersection& intersection,
|
|
const int insideElemIdx,
|
|
const int insideFaceIdx,
|
|
const int outsideElemIdx,
|
|
const int outsideFaceIdx,
|
|
DimVector& faceCenterInside,
|
|
DimVector& faceCenterOutside,
|
|
DimVector& faceAreaNormal,
|
|
/*isCpGrid=*/std::true_type) const
|
|
{
|
|
int faceIdx = intersection.id();
|
|
faceCenterInside = grid_.faceCenterEcl(insideElemIdx, insideFaceIdx);
|
|
faceCenterOutside = grid_.faceCenterEcl(outsideElemIdx, outsideFaceIdx);
|
|
faceAreaNormal = grid_.faceAreaNormalEcl(faceIdx);
|
|
}
|
|
|
|
template<class Grid, class GridView, class ElementMapper, class CartesianIndexMapper, class Scalar>
|
|
std::tuple<std::vector<NNCdata>, std::vector<NNCdata>>
|
|
EclTransmissibility<Grid,GridView,ElementMapper,CartesianIndexMapper,Scalar>::
|
|
applyNncToGridTrans_(const std::unordered_map<std::size_t,int>& cartesianToCompressed)
|
|
{
|
|
// First scale NNCs with EDITNNC.
|
|
std::vector<NNCdata> unprocessedNnc;
|
|
std::vector<NNCdata> processedNnc;
|
|
const auto& nnc_input = eclState_.getInputNNC().input();
|
|
if (nnc_input.empty())
|
|
return std::make_tuple(processedNnc, unprocessedNnc);
|
|
|
|
for (const auto& nncEntry : nnc_input) {
|
|
auto c1 = nncEntry.cell1;
|
|
auto c2 = nncEntry.cell2;
|
|
auto lowIt = cartesianToCompressed.find(c1);
|
|
auto highIt = cartesianToCompressed.find(c2);
|
|
int low = (lowIt == cartesianToCompressed.end())? -1 : lowIt->second;
|
|
int high = (highIt == cartesianToCompressed.end())? -1 : highIt->second;
|
|
|
|
if (low > high)
|
|
std::swap(low, high);
|
|
|
|
if (low == -1 && high == -1)
|
|
// Silently discard as it is not between active cells
|
|
continue;
|
|
|
|
if (low == -1 || high == -1) {
|
|
// Discard the NNC if it is between active cell and inactive cell
|
|
std::ostringstream sstr;
|
|
sstr << "NNC between active and inactive cells ("
|
|
<< low << " -> " << high << ") with globalcell is (" << c1 << "->" << c2 <<")";
|
|
OpmLog::warning(sstr.str());
|
|
continue;
|
|
}
|
|
|
|
auto candidate = trans_.find(isId(low, high));
|
|
|
|
if (candidate == trans_.end())
|
|
// This NNC is not resembled by the grid. Save it for later
|
|
// processing with local cell values
|
|
unprocessedNnc.push_back(nncEntry);
|
|
else {
|
|
// NNC is represented by the grid and might be a neighboring connection
|
|
// In this case the transmissibilty is added to the value already
|
|
// set or computed.
|
|
candidate->second += nncEntry.trans;
|
|
processedNnc.push_back(nncEntry);
|
|
}
|
|
}
|
|
return std::make_tuple(processedNnc, unprocessedNnc);
|
|
}
|
|
|
|
template<class Grid, class GridView, class ElementMapper, class CartesianIndexMapper, class Scalar>
|
|
void EclTransmissibility<Grid,GridView,ElementMapper,CartesianIndexMapper,Scalar>::
|
|
applyEditNncToGridTrans_(const std::unordered_map<std::size_t,int>& globalToLocal)
|
|
{
|
|
const auto& nnc_input = eclState_.getInputNNC();
|
|
const auto& editNnc = nnc_input.edit();
|
|
if (editNnc.empty())
|
|
return;
|
|
const auto& cartDims = cartMapper_.cartesianDimensions();
|
|
|
|
auto format_ijk = [&cartDims](std::size_t cell) -> std::string {
|
|
auto i = cell % cartDims[0]; cell /= cartDims[0];
|
|
auto j = cell % cartDims[1];
|
|
auto k = cell / cartDims[1];
|
|
|
|
return fmt::format("({},{},{})", i + 1,j + 1,k + 1);
|
|
};
|
|
|
|
auto print_warning = [&format_ijk, &nnc_input] (const NNCdata& nnc) {
|
|
const auto& location = nnc_input.edit_location( nnc );
|
|
auto warning = fmt::format("Problem with EDITNNC keyword\n"
|
|
"In {} line {} \n"
|
|
"No NNC defined for connection {} -> {}", location.filename,
|
|
location.lineno, format_ijk(nnc.cell1), format_ijk(nnc.cell2));
|
|
OpmLog::warning("EDITNNC", warning);
|
|
};
|
|
|
|
// editNnc is supposed to only reference non-neighboring connections and not
|
|
// neighboring connections. Use all entries for scaling if there is an NNC.
|
|
// variable nnc incremented in loop body.
|
|
auto nnc = editNnc.begin();
|
|
auto end = editNnc.end();
|
|
std::size_t warning_count = 0;
|
|
while (nnc != end) {
|
|
auto c1 = nnc->cell1;
|
|
auto c2 = nnc->cell2;
|
|
auto lowIt = globalToLocal.find(c1);
|
|
auto highIt = globalToLocal.find(c2);
|
|
|
|
if (lowIt == globalToLocal.end() || highIt == globalToLocal.end()) {
|
|
print_warning(*nnc);
|
|
++nnc;
|
|
warning_count++;
|
|
continue;
|
|
}
|
|
|
|
auto low = lowIt->second, high = highIt->second;
|
|
|
|
if (low > high)
|
|
std::swap(low, high);
|
|
|
|
auto candidate = trans_.find(isId(low, high));
|
|
if (candidate == trans_.end()) {
|
|
print_warning(*nnc);
|
|
++nnc;
|
|
warning_count++;
|
|
}
|
|
else {
|
|
// NNC exists
|
|
while (nnc!= end && c1==nnc->cell1 && c2==nnc->cell2) {
|
|
candidate->second *= nnc->trans;
|
|
++nnc;
|
|
}
|
|
}
|
|
}
|
|
|
|
if (warning_count > 0) {
|
|
auto warning = fmt::format("Problems with EDITNNC keyword\n"
|
|
"A total of {} connections not defined in grid", warning_count);
|
|
OpmLog::warning(warning);
|
|
}
|
|
}
|
|
|
|
template<class Grid, class GridView, class ElementMapper, class CartesianIndexMapper, class Scalar>
|
|
void EclTransmissibility<Grid,GridView,ElementMapper,CartesianIndexMapper,Scalar>::
|
|
computeHalfTrans_(Scalar& halfTrans,
|
|
const DimVector& areaNormal,
|
|
int faceIdx, // in the reference element that contains the intersection
|
|
const DimVector& distance,
|
|
const DimMatrix& perm) const
|
|
{
|
|
assert(faceIdx >= 0);
|
|
unsigned dimIdx = faceIdx/2;
|
|
assert(dimIdx < dimWorld);
|
|
halfTrans = perm[dimIdx][dimIdx];
|
|
|
|
Scalar val = 0;
|
|
for (unsigned i = 0; i < areaNormal.size(); ++i)
|
|
val += areaNormal[i]*distance[i];
|
|
|
|
halfTrans *= std::abs(val);
|
|
halfTrans /= distance.two_norm2();
|
|
}
|
|
|
|
template<class Grid, class GridView, class ElementMapper, class CartesianIndexMapper, class Scalar>
|
|
void EclTransmissibility<Grid,GridView,ElementMapper,CartesianIndexMapper,Scalar>::
|
|
computeHalfDiffusivity_(Scalar& halfDiff,
|
|
const DimVector& areaNormal,
|
|
const DimVector& distance,
|
|
const Scalar& poro) const
|
|
{
|
|
halfDiff = poro;
|
|
Scalar val = 0;
|
|
for (unsigned i = 0; i < areaNormal.size(); ++i)
|
|
val += areaNormal[i]*distance[i];
|
|
|
|
halfDiff *= std::abs(val);
|
|
halfDiff /= distance.two_norm2();
|
|
}
|
|
|
|
template<class Grid, class GridView, class ElementMapper, class CartesianIndexMapper, class Scalar>
|
|
typename EclTransmissibility<Grid,GridView,ElementMapper,CartesianIndexMapper,Scalar>::DimVector
|
|
EclTransmissibility<Grid,GridView,ElementMapper,CartesianIndexMapper,Scalar>::
|
|
distanceVector_(const DimVector& center,
|
|
int faceIdx, // in the reference element that contains the intersection
|
|
unsigned elemIdx,
|
|
const std::array<std::vector<DimVector>, dimWorld>& axisCentroids) const
|
|
{
|
|
assert(faceIdx >= 0);
|
|
unsigned dimIdx = faceIdx/2;
|
|
assert(dimIdx < dimWorld);
|
|
DimVector x = center;
|
|
x -= axisCentroids[dimIdx][elemIdx];
|
|
|
|
return x;
|
|
}
|
|
|
|
template<class Grid, class GridView, class ElementMapper, class CartesianIndexMapper, class Scalar>
|
|
void EclTransmissibility<Grid,GridView,ElementMapper,CartesianIndexMapper,Scalar>::
|
|
applyMultipliers_(Scalar& trans,
|
|
unsigned faceIdx,
|
|
unsigned cartElemIdx,
|
|
const TransMult& transMult) const
|
|
{
|
|
// apply multiplyer for the transmissibility of the face. (the
|
|
// face index is the index of the reference-element face which
|
|
// contains the intersection of interest.)
|
|
switch (faceIdx) {
|
|
case 0: // left
|
|
trans *= transMult.getMultiplier(cartElemIdx, FaceDir::XMinus);
|
|
break;
|
|
case 1: // right
|
|
trans *= transMult.getMultiplier(cartElemIdx, FaceDir::XPlus);
|
|
break;
|
|
|
|
case 2: // front
|
|
trans *= transMult.getMultiplier(cartElemIdx, FaceDir::YMinus);
|
|
break;
|
|
case 3: // back
|
|
trans *= transMult.getMultiplier(cartElemIdx, FaceDir::YPlus);
|
|
break;
|
|
|
|
case 4: // bottom
|
|
trans *= transMult.getMultiplier(cartElemIdx, FaceDir::ZMinus);
|
|
break;
|
|
case 5: // top
|
|
trans *= transMult.getMultiplier(cartElemIdx, FaceDir::ZPlus);
|
|
break;
|
|
}
|
|
}
|
|
|
|
template<class Grid, class GridView, class ElementMapper, class CartesianIndexMapper, class Scalar>
|
|
void EclTransmissibility<Grid,GridView,ElementMapper,CartesianIndexMapper,Scalar>::
|
|
applyNtg_(Scalar& trans,
|
|
unsigned faceIdx,
|
|
unsigned elemIdx,
|
|
const std::vector<double>& ntg) const
|
|
{
|
|
// apply multiplyer for the transmissibility of the face. (the
|
|
// face index is the index of the reference-element face which
|
|
// contains the intersection of interest.)
|
|
switch (faceIdx) {
|
|
case 0: // left
|
|
trans *= ntg[elemIdx];
|
|
break;
|
|
case 1: // right
|
|
trans *= ntg[elemIdx];
|
|
break;
|
|
|
|
case 2: // front
|
|
trans *= ntg[elemIdx];
|
|
break;
|
|
case 3: // back
|
|
trans *= ntg[elemIdx];
|
|
break;
|
|
|
|
// NTG does not apply to top and bottom faces
|
|
}
|
|
}
|
|
|
|
#ifdef HAVE_DUNE_FEM
|
|
template class EclTransmissibility<Dune::CpGrid,
|
|
Dune::GridView<Dune::Fem::GridPart2GridViewTraits<Dune::Fem::AdaptiveLeafGridPart<Dune::CpGrid, Dune::PartitionIteratorType(4), false>>>,
|
|
Dune::MultipleCodimMultipleGeomTypeMapper<Dune::GridView<Dune::Fem::GridPart2GridViewTraits<Dune::Fem::AdaptiveLeafGridPart<Dune::CpGrid, Dune::PartitionIteratorType(4), false>>>>,
|
|
Dune::CartesianIndexMapper<Dune::CpGrid>,
|
|
double>;
|
|
template class EclTransmissibility<Dune::CpGrid,
|
|
Dune::Fem::GridPart2GridViewImpl<
|
|
Dune::Fem::AdaptiveLeafGridPart<
|
|
Dune::CpGrid,
|
|
Dune::PartitionIteratorType(4),
|
|
false> >,
|
|
Dune::MultipleCodimMultipleGeomTypeMapper<
|
|
Dune::Fem::GridPart2GridViewImpl<
|
|
Dune::Fem::AdaptiveLeafGridPart<
|
|
Dune::CpGrid,
|
|
Dune::PartitionIteratorType(4),
|
|
false> > >,
|
|
Dune::CartesianIndexMapper<Dune::CpGrid>,
|
|
double>;
|
|
#if HAVE_DUNE_ALUGRID
|
|
#if HAVE_MPI
|
|
using ALUGrid3CN = Dune::ALUGrid<3, 3, Dune::cube, Dune::nonconforming, Dune::ALUGridMPIComm>;
|
|
#else
|
|
using ALUGrid3CN = Dune::ALUGrid<3, 3, Dune::cube, Dune::nonconforming, Dune::ALUGridNoComm>;
|
|
#endif //HAVE_MPI
|
|
|
|
template class EclTransmissibility<ALUGrid3CN, Dune::GridView<Dune::Fem::GridPart2GridViewTraits<Dune::Fem::AdaptiveLeafGridPart<ALUGrid3CN, Dune::PartitionIteratorType(4), false>>>,
|
|
Dune::MultipleCodimMultipleGeomTypeMapper<Dune::GridView<Dune::Fem::GridPart2GridViewTraits<Dune::Fem::AdaptiveLeafGridPart<ALUGrid3CN, Dune::PartitionIteratorType(4), false>>>>,
|
|
Dune::CartesianIndexMapper<ALUGrid3CN>,
|
|
double>;
|
|
template class EclTransmissibility<ALUGrid3CN,
|
|
Dune::Fem::GridPart2GridViewImpl<
|
|
Dune::Fem::AdaptiveLeafGridPart<
|
|
ALUGrid3CN,
|
|
Dune::PartitionIteratorType(4),
|
|
false> >,
|
|
Dune::MultipleCodimMultipleGeomTypeMapper<
|
|
Dune::Fem::GridPart2GridViewImpl<
|
|
Dune::Fem::AdaptiveLeafGridPart<
|
|
ALUGrid3CN,
|
|
Dune::PartitionIteratorType(4),
|
|
false> > >,
|
|
Dune::CartesianIndexMapper<ALUGrid3CN>,
|
|
double>;
|
|
#endif //HAVE_DUNE_ALUGRID
|
|
|
|
#else // !DUNE_FEM
|
|
template class EclTransmissibility<Dune::CpGrid,
|
|
Dune::GridView<Dune::DefaultLeafGridViewTraits<Dune::CpGrid>>,
|
|
Dune::MultipleCodimMultipleGeomTypeMapper<Dune::GridView<Dune::DefaultLeafGridViewTraits<Dune::CpGrid>>>,
|
|
Dune::CartesianIndexMapper<Dune::CpGrid>,
|
|
double>;
|
|
|
|
#if HAVE_DUNE_ALUGRID
|
|
#if HAVE_MPI
|
|
using ALUGrid3CN = Dune::ALUGrid<3, 3, Dune::cube, Dune::nonconforming, Dune::ALUGridMPIComm>;
|
|
#else
|
|
using ALUGrid3CN = Dune::ALUGrid<3, 3, Dune::cube, Dune::nonconforming, Dune::ALUGridNoComm>;
|
|
#endif //HAVE_MPI
|
|
|
|
template class
|
|
EclTransmissibility<ALUGrid3CN,
|
|
Dune::GridView<
|
|
Dune::ALU3dLeafGridViewTraits<const ALUGrid3CN, Dune::PartitionIteratorType(4)>
|
|
>,
|
|
Dune::MultipleCodimMultipleGeomTypeMapper<
|
|
Dune::GridView<Dune::ALU3dLeafGridViewTraits<const ALUGrid3CN, Dune::PartitionIteratorType(4)>>
|
|
>,
|
|
Dune::CartesianIndexMapper<ALUGrid3CN>,
|
|
double>;
|
|
#endif //HAVE_DUNE_ALUGRID
|
|
#endif //HAVE_DUNE_FEM
|
|
|
|
template class EclTransmissibility<Dune::PolyhedralGrid<3,3,double>,
|
|
Dune::GridView<Dune::PolyhedralGridViewTraits<3, 3, double, Dune::PartitionIteratorType(4)>>, Dune::MultipleCodimMultipleGeomTypeMapper<Dune::GridView<Dune::PolyhedralGridViewTraits<3,3,double,Dune::PartitionIteratorType(4)>>>,
|
|
Dune::CartesianIndexMapper<Dune::PolyhedralGrid<3,3,double>>,
|
|
double>;
|
|
|
|
} // namespace Opm
|