mirror of
https://github.com/OPM/opm-simulators.git
synced 2025-02-25 18:55:30 -06:00
move all applications into their top-level directory
thanks to [at]akva2 for the suggestion.
This commit is contained in:
Andreas Lauser
369
ebos/ecltransmissibility.hh
Normal file
369
ebos/ecltransmissibility.hh
Normal file
@@ -0,0 +1,369 @@
|
||||
// -*- 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.
|
||||
*/
|
||||
/*!
|
||||
* \file
|
||||
*
|
||||
* \copydoc Ewoms::EclTransmissibility
|
||||
*/
|
||||
#ifndef EWOMS_ECL_TRANSMISSIBILITY_HH
|
||||
#define EWOMS_ECL_TRANSMISSIBILITY_HH
|
||||
|
||||
#include <ewoms/common/propertysystem.hh>
|
||||
|
||||
#include <opm/parser/eclipse/EclipseState/EclipseState.hpp>
|
||||
#include <opm/parser/eclipse/EclipseState/Grid/GridProperties.hpp>
|
||||
#include <opm/parser/eclipse/EclipseState/Grid/FaceDir.hpp>
|
||||
#include <opm/parser/eclipse/EclipseState/Grid/TransMult.hpp>
|
||||
#include <opm/common/ErrorMacros.hpp>
|
||||
#include <opm/common/Exceptions.hpp>
|
||||
|
||||
#include <dune/common/version.hh>
|
||||
#include <dune/common/fvector.hh>
|
||||
#include <dune/common/fmatrix.hh>
|
||||
|
||||
#include <array>
|
||||
#include <vector>
|
||||
#include <unordered_map>
|
||||
|
||||
namespace Ewoms {
|
||||
namespace Properties {
|
||||
NEW_PROP_TAG(GridView);
|
||||
NEW_PROP_TAG(Scalar);
|
||||
NEW_PROP_TAG(Simulator);
|
||||
}
|
||||
|
||||
/*!
|
||||
* \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()
|
||||
{ update(); }
|
||||
|
||||
void update()
|
||||
{
|
||||
const auto& problem = simulator_.problem();
|
||||
const auto& gridManager = simulator_.gridManager();
|
||||
const auto& gridView = simulator_.gridView();
|
||||
const auto& elementMapper = simulator_.model().elementMapper();
|
||||
const auto& cartMapper = gridManager.cartesianIndexMapper();
|
||||
const auto eclState = gridManager.eclState();
|
||||
const auto eclGrid = eclState->getInputGrid();
|
||||
const auto& transMult = eclState->getTransMult();
|
||||
|
||||
const std::vector<double>& ntg =
|
||||
eclState->get3DProperties().getDoubleGridProperty("NTG").getData();
|
||||
|
||||
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
|
||||
|
||||
// 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.
|
||||
unsigned cartesianCellIdx = cartMapper.cartesianIndex(elemIdx);
|
||||
const auto& centroid = eclGrid.getCellCenter(cartesianCellIdx);
|
||||
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);
|
||||
|
||||
// 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 insideCartElemIdx = cartMapper.cartesianIndex(insideElemIdx);
|
||||
unsigned outsideCartElemIdx = cartMapper.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, insideCartElemIdx, ntg);
|
||||
applyNtg_(halfTrans2, outsideFaceIdx, outsideCartElemIdx, 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 ...
|
||||
applyMultipliers_(trans, insideFaceIdx, insideCartElemIdx, transMult);
|
||||
// ... and outside elements
|
||||
applyMultipliers_(trans, outsideFaceIdx, outsideCartElemIdx, transMult);
|
||||
|
||||
// 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 *= transMult.getRegionMultiplier(insideCartElemIdx,
|
||||
outsideCartElemIdx,
|
||||
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
|
||||
{
|
||||
Scalar isArea = is.geometry().volume();
|
||||
DimVector n = is.centerUnitOuterNormal();
|
||||
n *= isArea;
|
||||
|
||||
unsigned dimIdx = faceIdx/2;
|
||||
assert(dimIdx < dimWorld);
|
||||
halfTrans = perm[dimIdx][dimIdx];
|
||||
//halfTrans *= isArea;
|
||||
|
||||
Scalar val = 0;
|
||||
for (unsigned i = 0; i < n.size(); ++i)
|
||||
val += n[i]*distance[i];
|
||||
|
||||
halfTrans *= std::abs<Scalar>(val);
|
||||
halfTrans /= distance.two_norm2();
|
||||
}
|
||||
|
||||
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;
|
||||
}
|
||||
|
||||
void applyMultipliers_(Scalar& trans, unsigned faceIdx, unsigned cartElemIdx,
|
||||
const Opm::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, Opm::FaceDir::XMinus);
|
||||
break;
|
||||
case 1: // right
|
||||
trans *= transMult.getMultiplier(cartElemIdx, Opm::FaceDir::XPlus);
|
||||
break;
|
||||
|
||||
case 2: // front
|
||||
trans *= transMult.getMultiplier(cartElemIdx, Opm::FaceDir::YMinus);
|
||||
break;
|
||||
case 3: // back
|
||||
trans *= transMult.getMultiplier(cartElemIdx, Opm::FaceDir::YPlus);
|
||||
break;
|
||||
|
||||
case 4: // bottom
|
||||
trans *= transMult.getMultiplier(cartElemIdx, Opm::FaceDir::ZMinus);
|
||||
break;
|
||||
case 5: // top
|
||||
trans *= transMult.getMultiplier(cartElemIdx, Opm::FaceDir::ZPlus);
|
||||
break;
|
||||
}
|
||||
}
|
||||
|
||||
void applyNtg_(Scalar& trans, unsigned faceIdx, unsigned cartElemIdx,
|
||||
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[cartElemIdx];
|
||||
break;
|
||||
case 1: // right
|
||||
trans *= ntg[cartElemIdx];
|
||||
break;
|
||||
|
||||
case 2: // front
|
||||
trans *= ntg[cartElemIdx];
|
||||
break;
|
||||
case 3: // back
|
||||
trans *= ntg[cartElemIdx];
|
||||
break;
|
||||
|
||||
// NTG does not apply to top and bottom faces
|
||||
}
|
||||
}
|
||||
|
||||
const Simulator& simulator_;
|
||||
std::unordered_map<std::uint64_t, Scalar> trans_;
|
||||
};
|
||||
|
||||
} // namespace Ewoms
|
||||
|
||||
#endif
|
||||
Reference in New Issue
Block a user