mirror of
https://github.com/OPM/opm-simulators.git
synced 2025-02-25 18:55:30 -06:00
a63a4f2bdc
(instead of using 'int'.) This triggered quite a few compiler warnings which are also dealt-with by this patch.
381 lines
14 KiB
C++
381 lines
14 KiB
C++
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
|
|
// vi: set et ts=4 sw=4 sts=4:
|
|
/*
|
|
Copyright (C) 2014 by Andreas Lauser
|
|
|
|
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/>.
|
|
*/
|
|
/*!
|
|
* \file
|
|
*
|
|
* \copydoc Ewoms::EclTransmissibility
|
|
*/
|
|
#ifndef EWOMS_ECL_TRANSMISSIBILITY_HH
|
|
#define EWOMS_ECL_TRANSMISSIBILITY_HH
|
|
|
|
#include "eclgridmanager.hh"
|
|
|
|
#include <dune/common/version.hh>
|
|
#include <dune/common/fvector.hh>
|
|
#include <dune/common/fmatrix.hh>
|
|
|
|
#include <array>
|
|
#include <vector>
|
|
#include <unordered_map>
|
|
|
|
namespace Ewoms {
|
|
/*!
|
|
* \ingroup EclBlackOilSimulator
|
|
*
|
|
* \brief This class calculates the transmissibilites for grid faces according to the
|
|
* Eclipse Technical Description.
|
|
*/
|
|
template <class TypeTag>
|
|
class EclTransmissibility
|
|
{
|
|
typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
|
|
typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
|
|
typedef typename GET_PROP_TYPE(TypeTag, Simulator) Simulator;
|
|
typedef typename GridView::Intersection Intersection;
|
|
|
|
// Grid and world dimension
|
|
enum { dimWorld = GridView::dimensionworld };
|
|
|
|
typedef Dune::FieldMatrix<Scalar, dimWorld, dimWorld> DimMatrix;
|
|
typedef Dune::FieldVector<Scalar, dimWorld> DimVector;
|
|
|
|
public:
|
|
EclTransmissibility(const Simulator& simulator)
|
|
: simulator_(simulator)
|
|
{}
|
|
|
|
/*!
|
|
* \brief Actually compute the transmissibilty over a face as a pre-compute step.
|
|
*
|
|
* This code actually uses the direction specific "centroids" of
|
|
* each element. These "centroids" are _not_ the identical
|
|
* barycenter of the element, but the middle of the centers of the
|
|
* faces of the logical Cartesian cells, i.e., the centers of the
|
|
* faces of the reference elements. We do things this way because
|
|
* the barycenter of the element can be located outside of the
|
|
* element for sufficiently "ugly" (i.e., thin and "non-flat")
|
|
* elements which in turn leads to quite wrong
|
|
* permeabilities. This approach is probably not always correct
|
|
* either but at least it seems to be much better.
|
|
*/
|
|
void finishInit()
|
|
{
|
|
const auto& elementMapper = simulator_.model().elementMapper();
|
|
const auto& gridView = simulator_.gridView();
|
|
const auto& problem = simulator_.problem();
|
|
|
|
unsigned numElements = elementMapper.size();
|
|
|
|
// this code assumes that the DOFs are the elements. (i.e., an
|
|
// ECFV spatial discretization with TPFA). if you try to use
|
|
// it with something else, you're currently out of luck,
|
|
// sorry!
|
|
assert(simulator_.model().numGridDof() == numElements);
|
|
|
|
// 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>();
|
|
for (; elemIt != elemEndIt; ++elemIt) {
|
|
const auto& elem = *elemIt;
|
|
#if DUNE_VERSION_NEWER(DUNE_COMMON, 2,4)
|
|
unsigned elemIdx = elementMapper.index(elem);
|
|
#else
|
|
unsigned elemIdx = elementMapper.map(elem);
|
|
#endif
|
|
|
|
// get the geometry of the current element
|
|
const auto& geom = elem.geometry();
|
|
|
|
// compute the axis specific "centroids" used for the
|
|
// transmissibilities
|
|
for (unsigned dimIdx = 0; dimIdx < dimWorld; ++dimIdx) {
|
|
DimVector x0Local(0.5);
|
|
DimVector x1Local(0.5);
|
|
|
|
x0Local[dimIdx] = 0.0;
|
|
x1Local[dimIdx] = 1.0;
|
|
|
|
DimVector x = geom.global(x0Local);
|
|
x += geom.global(x1Local);
|
|
x /= 2;
|
|
|
|
axisCentroids[dimIdx][elemIdx] = x;
|
|
}
|
|
}
|
|
|
|
const auto& gridManager = simulator_.gridManager();
|
|
Opm::EclipseStateConstPtr eclState = gridManager.eclState();
|
|
auto multipliers = eclState->getTransMult();
|
|
const std::vector<double>& multx =
|
|
eclState->getDoubleGridProperty("MULTX")->getData();
|
|
const std::vector<double>& multy =
|
|
eclState->getDoubleGridProperty("MULTY")->getData();
|
|
const std::vector<double>& multz =
|
|
eclState->getDoubleGridProperty("MULTZ")->getData();
|
|
const std::vector<double>& multxMinus =
|
|
eclState->getDoubleGridProperty("MULTX-")->getData();
|
|
const std::vector<double>& multyMinus =
|
|
eclState->getDoubleGridProperty("MULTY-")->getData();
|
|
const std::vector<double>& multzMinus =
|
|
eclState->getDoubleGridProperty("MULTZ-")->getData();
|
|
|
|
const std::vector<double>& ntg =
|
|
eclState->getDoubleGridProperty("NTG")->getData();
|
|
|
|
// 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). unfortunately, this method is not available
|
|
// in GCC 4.4. (but I cannot say for sure when it became available so I play safe
|
|
// and limit it to GCC >= 4.8)
|
|
#if defined __clang__ || (__GNUC__ > 4 && __GNUC_MINOR__ >= 8)
|
|
trans_.reserve(numElements*3*1.05);
|
|
#endif
|
|
|
|
// compute the transmissibilities for all intersections
|
|
elemIt = gridView.template begin</*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;
|
|
|
|
// ignore boundary intersections for now (TODO?)
|
|
if (intersection.boundary())
|
|
continue;
|
|
|
|
const auto& inside = intersection.inside();
|
|
const auto& outside = intersection.outside();
|
|
#if DUNE_VERSION_NEWER(DUNE_COMMON, 2,4)
|
|
unsigned insideElemIdx = elementMapper.index(inside);
|
|
unsigned outsideElemIdx = elementMapper.index(outside);
|
|
#else
|
|
unsigned insideElemIdx = elementMapper.map(*inside);
|
|
unsigned outsideElemIdx = elementMapper.map(*outside);
|
|
#endif
|
|
|
|
// we only need to calculate a face's transmissibility
|
|
// once...
|
|
if (insideElemIdx > outsideElemIdx)
|
|
continue;
|
|
|
|
unsigned cartesianElemIdxInside = gridManager.cartesianIndex(insideElemIdx);
|
|
unsigned cartesianElemIdxOutside = gridManager.cartesianIndex(outsideElemIdx);
|
|
|
|
// local indices of the faces of the inside and
|
|
// outside elements which contain the intersection
|
|
unsigned insideFaceIdx = intersection.indexInInside();
|
|
unsigned outsideFaceIdx = intersection.indexInOutside();
|
|
|
|
Scalar halfTrans1;
|
|
Scalar halfTrans2;
|
|
|
|
computeHalfTrans_(halfTrans1,
|
|
intersection,
|
|
insideFaceIdx,
|
|
distanceVector_(intersection,
|
|
intersection.indexInInside(),
|
|
insideElemIdx,
|
|
axisCentroids),
|
|
problem.intrinsicPermeability(insideElemIdx));
|
|
computeHalfTrans_(halfTrans2,
|
|
intersection,
|
|
outsideFaceIdx,
|
|
distanceVector_(intersection,
|
|
intersection.indexInOutside(),
|
|
outsideElemIdx,
|
|
axisCentroids),
|
|
problem.intrinsicPermeability(outsideElemIdx));
|
|
|
|
applyNtg_(halfTrans1, insideFaceIdx, insideElemIdx, ntg);
|
|
applyNtg_(halfTrans2, outsideFaceIdx, outsideElemIdx, ntg);
|
|
|
|
// convert half transmissibilities to full face
|
|
// transmissibilities using the harmonic mean
|
|
Scalar trans = 1.0 / (1.0/halfTrans1 + 1.0/halfTrans2);
|
|
|
|
// apply the full face transmissibility multipliers
|
|
// for the inside ...
|
|
applyMultipliers_(trans, insideFaceIdx, insideElemIdx,
|
|
multx, multxMinus,
|
|
multy, multyMinus,
|
|
multz, multzMinus);
|
|
// ... and outside elements
|
|
applyMultipliers_(trans, outsideFaceIdx, outsideElemIdx,
|
|
multx, multxMinus,
|
|
multy, multyMinus,
|
|
multz, multzMinus);
|
|
|
|
// apply the region multipliers (cf. the MULTREGT keyword)
|
|
Opm::FaceDir::DirEnum faceDir;
|
|
switch (insideFaceIdx) {
|
|
case 0:
|
|
case 1:
|
|
faceDir = Opm::FaceDir::XPlus;
|
|
break;
|
|
|
|
case 2:
|
|
case 3:
|
|
faceDir = Opm::FaceDir::YPlus;
|
|
break;
|
|
|
|
case 4:
|
|
case 5:
|
|
faceDir = Opm::FaceDir::ZPlus;
|
|
break;
|
|
|
|
default:
|
|
OPM_THROW(std::logic_error, "Could not determine a face direction");
|
|
}
|
|
|
|
trans *= multipliers->getRegionMultiplier(cartesianElemIdxInside,
|
|
cartesianElemIdxOutside,
|
|
faceDir);
|
|
|
|
trans_[isId_(insideElemIdx, outsideElemIdx)] = trans;
|
|
}
|
|
}
|
|
}
|
|
|
|
Scalar transmissibility(unsigned elemIdx1, unsigned elemIdx2) const
|
|
{ return trans_.at(isId_(elemIdx1, elemIdx2)); }
|
|
|
|
private:
|
|
std::uint64_t isId_(unsigned elemIdx1, unsigned elemIdx2) const
|
|
{
|
|
static const unsigned elemIdxShift = 32; // bits
|
|
|
|
unsigned elemAIdx = std::min(elemIdx1, elemIdx2);
|
|
std::uint64_t elemBIdx = std::max(elemIdx1, elemIdx2);
|
|
|
|
return (elemBIdx<<elemIdxShift) + elemAIdx;
|
|
}
|
|
|
|
void computeHalfTrans_(Scalar& halfTrans,
|
|
const Intersection& is,
|
|
unsigned faceIdx, // in the reference element that contains the intersection
|
|
const DimVector& distance,
|
|
const DimMatrix& perm) const
|
|
{
|
|
unsigned dimIdx = faceIdx/2;
|
|
assert(dimIdx < dimWorld);
|
|
halfTrans = perm[dimIdx][dimIdx];
|
|
halfTrans *= is.geometry().volume();
|
|
|
|
const auto &normal = is.centerUnitOuterNormal();
|
|
Scalar val = 0;
|
|
for (unsigned i = 0; i < normal.size(); ++i)
|
|
val += is.centerUnitOuterNormal()[i]*distance[i];
|
|
|
|
halfTrans *= std::abs<Scalar>(val);
|
|
halfTrans /= distance*distance;
|
|
}
|
|
|
|
DimVector distanceVector_(const Intersection& is,
|
|
unsigned faceIdx, // in the reference element that contains the intersection
|
|
unsigned elemIdx,
|
|
const std::array<std::vector<DimVector>, dimWorld>& axisCentroids) const
|
|
{
|
|
unsigned dimIdx = faceIdx/2;
|
|
assert(dimIdx < dimWorld);
|
|
DimVector x = is.geometry().center();
|
|
x -= axisCentroids[dimIdx][elemIdx];
|
|
|
|
return x;
|
|
}
|
|
|
|
template <class MultScalar>
|
|
void applyMultipliers_(Scalar &trans, unsigned faceIdx, unsigned elemIdx,
|
|
const std::vector<MultScalar>& multx,
|
|
const std::vector<MultScalar>& multxMinus,
|
|
const std::vector<MultScalar>& multy,
|
|
const std::vector<MultScalar>& multyMinus,
|
|
const std::vector<MultScalar>& multz,
|
|
const std::vector<MultScalar>& multzMinus) 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 *= multxMinus[elemIdx];
|
|
break;
|
|
case 1: // right
|
|
trans *= multx[elemIdx];
|
|
break;
|
|
|
|
case 2: // front
|
|
trans *= multyMinus[elemIdx];
|
|
break;
|
|
case 3: // back
|
|
trans *= multy[elemIdx];
|
|
break;
|
|
|
|
case 4: // bottom
|
|
trans *= multzMinus[elemIdx];
|
|
break;
|
|
case 5: // top
|
|
trans *= multz[elemIdx];
|
|
break;
|
|
}
|
|
}
|
|
|
|
template <class NtgScalar>
|
|
void applyNtg_(Scalar &trans, unsigned faceIdx, unsigned elemIdx,
|
|
const std::vector<NtgScalar>& 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
|
|
}
|
|
}
|
|
|
|
const Simulator& simulator_;
|
|
std::unordered_map<std::uint64_t, Scalar> trans_;
|
|
};
|
|
|
|
} // namespace Ewoms
|
|
|
|
#endif
|